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SECTION  1 
INTRODUCTION 


An  overview  of  the  Doppler  satellite  tracking  data  processing  system  with  which 
this  document  is  largely  concerned  is  shown  in  Figure  1-1.  In  this  figure  the 
system's  data  processing  steps  are  represented  by  circles  and  locations  which 
provide  support  data  or  use  the  processed  results  are  represented  by  rectangles. 
Also  shown  is  the  basic  data  flow  between  processing  steps.  Note  that  only 
those  processing  steps  enclosed  by  the  dashed  box  are  performed  by  the  PULSAR 
data  editor. 


The  following  sections  discuss  in  detail  both  the  mathematical  and  the  pro¬ 
cessing  methodologies  employed  by  those  processing  steps  enclosed  within  the 
PULSAR  functional  boundary.  Included  in  this  discussion  are  descriptions  of 
pertinent  coordinate  reference  systems;  the  force  model  and  variational  eq¬ 
uations  used  by  PULSAR;  the  time  calibration  procedures;  observation  correction 
methods;  the  state  estimation  and  parameter  improvement  process;  data  quality 
statistics;  as  well  as  numerical  techniques  applied  during  processing.  Since 
the  analysis  reporting  function  is  essentially  a  bookkeeping  procedure  and  is 
non-mathematical  in  nature,  it  is  not  discussed  within  this  document.  It  should 
be  mentioned  that  the  PULSAR  program  is  an  outgrowth  of  the  CELEST^  program  used 
at  NSWC  for  many  years  and,  as  a  result,  both  programs  share  common  fundamental 
algorithms.  Although  this  document  is  self-contained,  the  reader  may  wish  to 
consult  reference  1  for  supporting  background  material. 
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SECTION  2 

COORDINATE  REFERENCE  SYSTEMS  AND  TRANSFORMATIONS 


A  number  of  coordinate  reference  frames  are  used  by  the  PULSAR  data  editing 
system.  These  include  the  following: 


(i)  the  basic  inertial  system; 

(ii)  the  mean  of  date  system; 

(iii)  the  true  of  date  system; 

( i v )  the  instantaneous  earth-fixed  system; 

(v)  the  earth-fixed  reference  ellipsoid  system; 

(vi)  the  geographical  coordinate  system;  and 

(vii)  the  time  of  closest  approach  (tca)  reference  system. 


Each  of  these  reference  systems,  as  well  as  transformations  between  them, 
will  be  discussed  in  some  detail  in  the  following  subsections. 


2 . 1  EQUATORIAL  INERTIAL  COORDINATE  SYSTEMS 

Reference  systems  (i),  (ii),  and  (iii)  above  are  categorized  as  equatorial 
inertial  systems.  However,  before  undertaking  a  description  of  those  inertial 
coordinate  systems,  it  is  necessary  to  define  certain  terms  of  reference. 
Fundamental  to  understanding  those  coordinate  frames  is  the  concept  of  the 
celestial  sphere.  This  is  an  imaginary  sphere  of  infinite  radius  the  "center" 
of  which  is  coincident  with  the  center  of  the  earth.  The  basic  components 
of  a  system  ■  . i ng  the  celestial  sphere  are  the  celestial  equator, 
celestial  pok  ’iptic  and  equinox.  The  celestial  equator  and  celestial 
poles  are  obtainc  majecting  the  earth's  equator  and  poles  onto  the  cele¬ 

stial  sphere.  The  projection  of  the  apparent  path  of  the  sun  onto  the  celestial 
sphere  forms  a  great  circle  inclined  to  the  celestial  equator  at  an  angl^  of 
approximately  23.5°.  This  great  circle  is  called  the  ecliptic  and  the  incli¬ 
nation  angle  is  called  the  obliquity  of  the  ecliptic.  The  ecliptic  intersects 
the  celestial  equator  at  two  points  called  the  vernal  equinox  and  the  autumnal 
equinox,  i.e.  the  ascending  and  descending  nodes,  respectively,  of  the  sun's 
apparent  orbit.  For  historical  reasons  the  vernal  equinox  is  commonly  refer¬ 
red  to  as  the  first  point  of  Aries  and  is  labeled  with  the  symbol  T.  A 
depiction  of  the  celestial  sphere  and  it's  components  is  given  in  Figure  2-1. 
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Table  2-2.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (B1950.0) 


-i 

a  ■ 

A. 

b. 

-c 

d  ■ 

A 

1. 

12.11684352 

-0.0529538653 

0. 15595 (-11) 

0 . 4600 (-19) 

200.01154773 

1.9712946873 

0.45354(-12) 

0 . 2800 (-19) 

24.23368705 

-0.1059077306 

0.31 190 (-11) 

0 . 9200 (-19) 

126.73383654 

26.3527929915 

-0. 16941 (-11) 

0. 1 060 ( - 18 ) 

357.92523074 

0.9856002627 

-0.1 1576 ( -12 ) 

-0. 6800( -19) 

214.53145564 

13.0649926982 

0.69061 ( -11 ) 

0 . 2950 ( -18) 

7. 

197.93677848 

2.9568949501 

0.33776(-12) 

-0.4000( -19) 

8. 

114.61699301 

26.4057468568 

-0.32537 ( -11 ) 

0. 6000( -19) 

9. 

341.26529219 

39.4177856897 

0. 52119( -11 ) 

0 . 4010 ( -18 ) 

10. 

202.08631698 

0.9856944245 

0. 56929( -12) 

0 . 9600 (-19) 

11. 

287.80916684 

-11.3165056059 

0 . 90538 (-11) 

0. 2170( -18) 

12. 

187.89470420 

2.0242485526 

-0.1 1059 (-11) 

-0. 1800( -19) 

13. 

272.20238089 

13.2878002932 

-0.86003(-ll) 

-0. 1890 ( - 18 ) 

14. 

286.72228880 

24.3814983042 

-0 . 21477 (-11) 

0. 7800( -19) 

15. 

266.64829917 

13.0120388329 

0.84656( -11 ) 

0 . 3410 ( -18) 

16. 

157.58538787 

-13.1179465635 

-0. 53466 (-11) 

-0 . 2490 ( -18) 

17. 

198.92466969 

37.6692985975 

-0. 1 0748 (-10) 

-0. 1 1 10 ( - 18 ) 

18. 

45.55408171 

0.2757614603 

-0. 1 7  066 (-10) 

-0. 5300 (-18) 

19. 

142.34062249 

1.7484870922 

0. 1 5960 ( -10) 

0. 5120( -18) 

20. 

329.14844866 

39.4707395551 

0 . 36524 (-11) 

0. 3550( -18) 

21. 

53.45612534 

50.7342912957 

-0 . 38419 ( -11) 

0. 1840 ( - 18 ) 

22. 

69.06291129 

26.1299853965 

0. 1 3812 (-10) 

0 . 5900 ( - 18 ) 

23. 

54.54300338 

15.0362873855 

0 . 73596 (-11) 

0 . 3230 ( -18) 

24. 

195.79674734 

52.4827783880 

0. 121 18( -10) 

0. 6960( -18) 

25. 

102.50014948 

26.4587007221 

-0 . 48132 (-11) 

0 . 1400 ( 

The  quantities  in 

the  parentheses  are 

the  exponents  of  10. 
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Table  2-1.  Coefficients  of  the  Nutation  Sine  and  Cosine  Series  (B1950.0) 
(Con ' t) 


:x  •  +  3  •  T 

i  t 

Y-  +  6  - T 
-C  A. 

-0. UK-06) 

0.0 

0.11K-06) 

0.0 

. 

0.11K-06) 

0.0 

-0 . Ill ( -06) 

0.55(-07) 

. 

-0.83(-07) 

0.55(-07) 

. 

0.83(-07 ) 

-0.551-07) 

-0 . 83 ( -07 ) 

0.0 

0.83(-07) 

0.0 

. 

-0.83(-07) 

0.0 

. 

-0 . 83 ( -07 ) 

0.0 

-0.55(-07) 

0.0 

. 

-0.551-07) 

0.0 

. 

-0.551-07) 

0.0 

0.551-07) 

0.0 

. 

-0.551-07) 

0.0 

.  I  -0 . SC ( -07 ) 

0.0 

!  -  V  51-07) 

0.0 

•  •' -07) 

0.0 

-  0 . 5  n  :  -  C  7  ) 

0.0 

0.0 

0.0 

! 

0.0 

i 

i 

0.0 

i 

! 

0.0 

0.0 

j 

0.0 

‘■ff  t  ionts  numbered  70  through  106  are  all  zero  for  this 
!  ::  O. •!)  version  of  the  transformation. 
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Table  2-1.  Coefficients  of  the  Nutation  Sine  and  Cosine  Series  (B1950.0) 
(con't) 


-i 

\ 

\  +  v 

26. 

-0.583(-06) 

0.0 

27. 

0.527 ( -06) 

-0.277 ( -06) 

28. 

0.443(-06)  -  0.3(-08)T 

0.0 

29. 

-0 . 416 ( -06 ) 

0.222( -06) 

30. 

-0.415(-06)  +  0.3(-08)T 

0.194(-06) 

31. 

0. 388( -06) 

-0. 194( -06) 

32. 

-0.361(-06) 

0. 194(-06) 

33. 

0.277(-06) 

0.0 

34. 

-0.277 ( -06) 

0. 138( -06) 

35. 

-0 . 250 ( -06 ) 

0. 138(-06) 

36.  I 

-0. 194 ( -06 ) 

0.0 

37. 

0.1 94 ( -06 ) 

-0.83( -07 ) 

38. 

0.1 66 ( -06 ) 

0.0 

39. 

-0 . 166 ( -06 ) 

0.83(-07) 

40. 

-0.1 66 ( -06 ) 

0.83(-07) 

41. 

-0. 166( -06) 

0 . 83 ( -07 ) 

42. 

0. 166( -06) 

-0. 55(-07 ) 

43. 

-0. 138( -06) 

0.83(-07) 

44. 

-0. 138( -06) 

0.83(-07) 

45. 

-0. 1 38 ( -06 ) 

0.83(-07) 

46. 

-0. 138( -06) 

0.83(-07 ) 

47. 

0. 1 38 ( -06 ) 

-0.83( -07 ) 

48. 

-0. UK-06) 

0.55(-07) 

49. 

0.11K-06) 

-0.55( -07) 

50. 

-0. UK-06) 

0.0 
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Table  2-1.  Coefficients  of  the  Nutation  Sine  and  Cosine  Series  (B1950.0) 


i 

i— 

+ 

v  *  siT 

1. 

-0.4789274(~02)  -0.4825( -05)T 

+0.2558460(-02)  +  0.252(-06)T 

2. 

-0.353601(-03)  -0 . 36 ( -07 )  T 

0. 153348( -03)  -  0.80(-07)T 

3. 

0.58002 ( -04)  +0. 5( -08)T 

-0 . 251 05 ( -04 )  +  0,11 ( -07)T 

4. 

-0. 56586(-04)  -0.5(-08)T 

0. 24548( -04)  -  0.14(-07)T 

5. 

0. 34984(-04)  -0.86(-07)T 

0.0 

6. 

0. 18751 ( -04)  +0.27(-08)T 

0.0 

7. 

-0. 13788 ( -04 )  +0.33(-07)T 

0.5991 ( -05)  -  0 . 17 ( -07 )T 

8. 

-0.9505(-05)  -0.11(-07)T 

0. 5083( -05) 

9. 

-0 . 7250 ( -05 ) 

0 . 3137 ( -05 )  -  0.27(-08)T 

10. 

0. 5937 ( -05)  -0. 1 3 ( -07 ) T 

-0 . 2  579 ( -05 )  +  0.83( -08)T 

11. 

-0.41 38 ( -05 ) 

0.0 

12. 

0 . 3445( -05)  +0.32(-08)T 

-0 . 1833 ( -05 ) 

13. 

0 . 3166 ( -05 ) 

-0 . 1388 ( -05) 

14. 

0. 1 666 ( -05 ) 

0.0 

15. 

0. 1611 ( -05) 

-0.861 (-06) 

16. 

-0. 1 583 ( -05) 

0.833(-06) 

17. 

-0.1 144 (-05) 

0. 61 1 ( -06) 

18. 

0. 1 2  50 ( -05) 

-0. 666( -06) 

19. 

0. 12 50 (-05) 

0.0 

20. 

-0. 1222 ( -05 ) 

0.638(-06) 

2i . 

-0.088 (-06) 

0. 388( -06) 

22. 

0.777 ( -06) 

0.0 

23. 

0.722 (-06) 

-0.305(-06) 

24. 

-  0 .722  1-06) 

0.305(-06) 

25 

••'(- 36) 

0.0 

Mi)  -  2,433,282.432,459,1)  /  36,525.0 
!  fie  quantities  in  the  parentheses  are  the  exponents  of  10 
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2.1.4  THE  TRUE  OF  DATE  INERTIAL  REFERENCE  SYSTEM  AND  THE  C  TRANSFORMATION  (con't) 


eC>  - 


23.44578773  -  0.356304  x  10_6(t-to)  -  0.66319  x  10'15(t-tQ)2 
+  0.10318  x  10‘19  (t-tQ)3  ,  B1950.0  _or_ 


23.43929111  -  0.3560346794  x  10  (t-t  )  -  0.122848274 

o 

x  10"15(t-to)2  +  0.103353  x  10-19(t-to)3  ,  J2000.0 


A^( ° )  = 


106 

Z  K  *  6iT> si"  U +  -  V  +  c^(t  -  lo>2  +di 

i=l  * 


1 0r  ,  o) 


lUb 

Ac(°)  =  £  hi  +  <5J)  cos  +  bjt  -  tQ)  +  c.(t  -  tQ)2  +  d.(t  -  to)3J. 

Equations  (9)  -  (10)  apply  for  both  the  B1950.0  and  J2000.0  epochs  and 


t  -  t0 
36525.' 0 


Values  for  the  a;,  B-,  y  ■,  and  6.  coefficients  are  given  in  Tables  2-1  and 
2-3.  Values  for  the  a^,  b^,  c^,  and  d^  coefficients  are  given  in  Tables 
2-2  and  2-4. 


For  an  arbitrary  vector  x»  one  has  the  following  relationships 


X"  =  C  X  '  =  C  D  x  , 
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e  =  TRUE  OBLIQUITY  OF  DATE 

i  =  MEAN  OBLIQUITY  OF  DATE 

Ae  =  NUTATION  IN  OBLIQUITY 

±4  =  NUTATION  IN  CELESTIAL  LONGITUDE 

AH  =  TAN*1  ITAN  A 0  COSe] 


FIGURE  2-4.  Nutation  Geometry  Showing  A e  and  Corrections 
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2.1.4  THE  TRUE  OF  DATE  INERTIAL  REFERENCE  SYSTEM  AND  THE  C  TRANSFORMATION 
The  true  of  date  system  is  the  inertial  reference  system  at  the  epoch  of  date 
which  has  its  (X,  Y,  Z)  axes  defined  as  follows: 

An 

X  is  the  unit  vector  along  the  true  vernal  equinox  of  date; 

A  it 

Z  is  the  unit  vector  along  the  earth's  true  axis  of  rotation 
of  date  and  is  positive  in  the  northern  hemisphere;  and 

Ah 

Y  is  the  unit  vector  which  completes  the  right-handed  orthogonol 

A"  An 

system  with  X  and  Z  . 

The  shift  from  the  mean  of  date  system  (as  defined  by  precessional  motion  only) 
to  the  true  of  date  system  is  resolved  into  nutation  induced  corrections  known 
as  the  nutation  in  ecliptic  longitude  (A\i>)  and  the  nutation  in  obliquity  (Ac). 
The  nutation  geometry  of  Figure  2-4  defines  these  corrections. 


The  transformation  from  the  mean  of  date  system  to  the  true  of  date  system  is 
accomplished  through  application  of  the  C  matrix  given  by: 


C  = 


^  cosAip  -sinA^cose  -sinA4>sine 

cosesinAijj  cosccosA^cose+sinesine  cosecosA^sine-sinecose 
\  sinesinAijj  sinecosA^cose  -  coscsine  sinecosA^sine+  cosecose 


/ 


,  (6) 


where 


e  -  e  +  Ac 


(7) 
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2.1.3  THE  MEAN  OF  DATE  INERTIAL  REFERENCE  SYSTEM 
'  0.15242971  x  10'4  (t-t  )  -  0.88704407  x 
-  0.23849  x  10'18  (t-tQ)3,  B1950.0 


«-( °)  =  S 


AND  THE  D  TRANSFORMATION  (con1 
10"13  (t-tQ)2 
or  , 


0.15243067  x  10‘4  (t-tQ)  -  0.88835960  x  10-13  (t-tQ)2 
•  -  0.23847663  x  10'18  (t-tQ)3,  J2000.0 


t) 


In  the  above  expressions 


t  = 


and 


Julian  Date  (JD),  B1950.0 

Julian  Ephemeris  Date  (JED),  J2000.0 


t 


o 


2433282.4234591,  B1950.0 

2451545.0  ,  J2000.0 


Thus  one  may  write  for  an  arbitrary  vector  ^  : 


(3) 


X*  =  ?  X  (4) 

where  x'  is  the  vector  in  the  mean  of  date  system  and  x  is  the  vector  in  the 
basic  reference  system.  Note  that  due  to  the  orthogonality  of  the  D  matrix 
one  has 

X  *  ?TJ'  .  (5) 
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FIGURE  2-3.  Precession  Geometry  Showing  j  ,  z  and  e 


TR  82-391 


2.1.3  THE  MEAN  OF  DATE  INERTIAL  REFERENCE  SYSTEM  AND  THE  D  TRANSFORMATION  (con't) 

The  difference  in  orientation  of  axes  of  the  basic  and  of  date  inertial 
systems  is  a  result  of  the  general  precession  which  is  discussed  in  sub¬ 
section  2.1.1.  The  complex  motion  of  the  general  precession  can  be  specified 
by  the  three  angles  £,  Z,  and  £.  These  angles  are  defined  by  the  geometry 
shown  in  Figure  2-3. 


The  transformation  from  the  mean  equator  and  equinox  of  either  of  the 
standard  epochs  (B1950.0  or  J2000.0)  to  the  mean  equator  and  equinox  of 
date  is  given  by 


D  = 


cosZcosQcos^-sinZsins 

sinZcosOcosc+cosZsin^ 

sindcoss 


-  cosZcosQsinC-sinZcosC 

-  sinZcos0sin£+cosZcos£ 

-  sin&sinc 


-  cosZsin& 

-  sinZsinQ 

cos& 


\ 

/ 


.  u) 


where 


/ 


< 


V 


f 


*C)  -< 


0.17529829  x  10  4  (t-tQ)  +  0.62884345  x  10'13(t-to)2 
+  0.10261  X  10-18(t-to)3,  B1950.0  or 

0.17539114  x  10-4  (t-tQ)  +  0.62856673  x  10"13(t-to)2 

+  0.10260087  x  10-18(t-to)3,  J2000.0 

0.17529829  x  10'4  (t-tQ)  +  0.22759135  x  10'12(t-to)2 

lO  O 

+  0.10261  x  10-1°  (t-t0)J,  B1950.0  .or. 

0.17539114  x  10'4  (t-tQ)  +  0.22793143  x  10'12(t-tQ)2 
+  0.10376951 x  10'18(t-to)3,  02000.0 


(2) 

(con’t  next 
page) 


1  Explanatory  Supplement  to  the  Astronomical  Ephemeris  and  the  American 
Ephemeris  and  Nautical  Almanac,  Her  Majesty's  Stationery  Office,  1961. 
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2.1.2  PULSAR  BASIC  INERTIAL  REFERENCE  SYSTEM 

When  integrating  the  satellite  equations  of  motion  (see  SECTION  5),  it  is 
necessary  to  select  a  reference  frame  which  is  free  of  all  movement.  Per 
the  above  discussion,  this  implies  that  the  integration  space  must  be  an 
inertial  system  fixed  to  a  given  equator  and  equinox,  i.e.  a  standard  epoch. 
PULSAR  may  optionally  adopt  as  its  standard  epoch  either  the  Besselian 
epoch  B1950.0  or  the  Julian  epoch  J2000.0  so  that  the  basic  inertial  refer¬ 
ence  frame  is  that  set  of  inertial  rectangular  coordinates  associated  with 
either  the  mean  equator  and  equinox  of  B1950.0  or  J2000.0,  respectively. 
Thus  the  (X,  Y,  Z)  axes  of  the  basic  inertial  reference  system  are  defined 
as  follows: 

A 

X  is  the  unit  vector  pointing  toward  the  mean  vernal  equinox 
of  B1950.0  or  J2000.0; 

A 

Z  is  the  unit  vector  pointing  along  the  earth's  mean  axis  of 
rotation  of  B1950.0  or  J2000.0  and  is  positive  in  the  northern 
hemisphere;  and 
A 

Y  is  the  unit  vector  which  completes  the  right-handed  orthogonal 

A  A 
system  with  X  and  Z. 

2.1.3  THE  MEAN  OF  DATE  INERTIAL  REFERENCE  SYSTEM  AND  THE  D  TRANSFORMATION 
The  mean  of  date  inertial  reference  system  differs  from  the  basic  inertial 
reference  system  only  in  that  its  reference  epoch  can  be  any  time,  e.g.  an 
observation  time.  Therefore  the  (X',  Y ' ,  •£' )  axes  of  this  reference  frame 
are  defined  as  follows: 

A 

X'  is  the  unit  vector  pointing  toward  the  mean  vernal  equinox  of  date; 

A 

Z1  is  the  unit  vector  pointing  along  the  earth's  mean  axis  of  rotation 
of  date  and  is  positive  in  the  northern  hemisphere;  and 

A 

Y'  is  the  unit  vector  which  completes  the  right-handed  orthogonal 
A  A 

system  with  X '  and  V  . 


L 
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FIGURE  2-2.  P recessional  and  Mutational  Motion 
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2.1  EQUATORIAL  INERTIAL  COORDINATE  SYSTEMS  (con't) 

All  equatorial  inertial  coordinate  systems  are  formulated  using  the  celes¬ 
tial  sphere  and  its  components  and  are  based  u,.cn  the  set  of  mutually 
orthogonal  axes  (X,  Y,  Z),  where  X  is  an  axis  aligned  in  the  direction  of 
the  vernal  equinox,  Z  is  aligned  in  the  direction  of  the  celestial 

pole,  and  Y  is  aligned  in  a  manner  to  complete  the  right-handed  ortho¬ 
gonal  system.  Unfortunately,  this  system  of  coordinates  is  not  fixed  in 
time,  but  is  time  varying.  This  is  due  to  the  fact  that  the  equinox  under¬ 
goes  a  slow  secular  movement  known  as  precession  and  a  small  periodic  varia¬ 
tion  known  as  nutation.  The  complications  associated  with  this  motion  can 
be  compensated  for.  This  is  the  subject  of  the  following  subsections. 

2.1.1  PRECESSION  AND  NUTATION 

The  gravitational  attractions  of  the  moon  and  sun  upon  the  earth's  equatorial 
bulge  tend  to  pull  the  earth's  equatorial  plane  into  coincidence  with  the 
plane  of  the  ecliptic.  However,  due  to  the  earth's  rotation,  these  attractions 
result  instead  in  a  conical  motion  of  the  earth's  axis  about  the  pole  of  the 
ecliptic  plane  (i.e.  the  intersection  with  the  celestial  sphere  of  an  imaginary 
line  projected  outwards  from  the  center  of  the  sphere  in  a  direction  perpen¬ 
dicular  to  the  plane  of  the  ecliptic).  This  motion  is  in  the  opposite  direction 
to  the  earth's  rotation  and  has  a  period  of  approximately  26,000  years. 

This  conical  motion  is  usually  treated  as  the  sum  of  two  components.  One 
consists  of  the  uniform  westward  motion  of  the  vernal  equinox  along  the  ecliptic. 
This  is  called  luni-solar  precession.  The  other  component  of  motion  is  per¬ 
iodic  and  is  called  nutation.  It  describes  the  departure  of  the  actual  position 
of  the  pole  of  the  equator  from  its  mean  position  as  defined  by  precessional 
motion  only. 

There  are  also  smaller  precessional  motions  of  the  ecliptic  plane  caused  by 
the  planetary  gravitational  attraction  and  general  relativistic  effects. 

These  motions  are  called  planetary  precession  and  geodesic  precession,  resp¬ 
ectively.  The  combined  effect  of  luni-solar  precession, planetary  precession 
and  geodesic  precession  is  called  general  precession.  The  superposition  of 
precessional  and  nutational  motion  is  depicted  in  Figure  2-2. 
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Table  2-2.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (B1950.0) 
(con’t) 


i. 

a . 

4. 

b. 

4. 

c  - 

4. 

d. 

4, 

26. 

175.77786068 

2.0772024179 

-0.26655(-ll) 

-0. 6400(-19) 

27. 

260.08553736 

13.3407541586 

-0. 10159(-10) 

-0.2350(-18) 

28. 

355.85046149 

1.9712005255 

-0.23153(-12) 

-0. 1360(-18) 

29. 

10.04207427  j 

0.9326463974 

0 . 14437 (-11) 

-0.2200(-19) 

30. 

195.86200923 

3.9424952128 

0 . 22200 ( -12 ) 

-0 . 1080 ( -18 ) 

31. 

84.30767668 

11.2635517406 

-0 . 74943 (-11) 

-0.1710(-18) 

32. 

153.37058798 

37.3935371371 

0.63179(-11) 

0.4190(-18) 

33. 

326.56276181 

-  0.3287153256 

0. 18625 ( -10) 

0. 5760(-18) 

34. 

14.19161277 

-  1.0385541281 

0. 16752 ( -11 ) 

0. 1140( -18) 

35. 

260.08553736 

13.3407541586 

-0. 10159(-10) 

-0. 2350( -18) 

36. 

285.73439759 

-10.330905343 

0 . 89381 ( -11 ) 

0.1490(-18) 

37. 

124.65906729 

27.3383932543 

-0.18099(-11) 

0.3800(-19) 

38. 

355.78520010 

50.5114837007 

0.11664(-10) 

0.6680(-18) 

39. 

298.83913233 

24.3285444389 

-0.58820(-12) 

0. 1240(-18) 

40. 

128.80860579 

25.3671927287 

-0. 15784 (-11) 

0. 1740( -18) 

41. 

267.98758099 

63.7992839940 

0 . 30642 (-11) 

0.4790(-18) 

42. 

269.07445903 

28.1012800838 

0. 14265 ( -10) 

0. 6180(-18) 

43. 

229.77622103 

-  1.8014409576 

-0. 14400(-10) 

-0.4660(-18) 

44. 

189.96947346 

1.0386482898 

-0. 99022( -12) 

0 . 5000( -19 ) 

45. 

85.39455472 

-24.4344521695 

0. 37072(-ll) 

-0.3200(-19) 

46. 

41.33928181 

50.7872451610 

-0 . 54014 ( -11) 

0 . 1 380 ( - 18 ) 

47. 

42.42615985 

15.0892412508 

0 . 58001 (-11) 

0 . 2770( -18) 

48. 

192.04424271 

0.0530480270 

-0.87445( -12) 

0. 1 180 ( - 18 ) 

49. 

154.45746602 

1.6955332269 

0. 17519( -10) 

0 . 5580 ( -18) 

50. 

143.36114440 

12.1907491521 

-0.1 07  38 (-11) 

0 . 3900( -19) 
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Table  2-2.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (B1950.0) 
(Con't) 


i 

a. . 

b. 

c  - 

d. 

A, 

■L 

A. 

A, 

51. 

71.20294194 

-23.3958980414 

0 . 20319 (-11) 

-0. 1460(-18) 

52. 

216.60622490 

12.0793924354 

0.70219(-11) 

0 . 3630 ( - 18 ) 

53. 

112.03130616 

-13.3937080239 

0 . 117 19 ( -10 ) 

0.2810(-18) 

54. 

183.67990431 

52.5357322533 

0. 10558(-10) 

0.6500(-18) 

55. 

57.67092524 

0.2228075950 

1 

-0. 15506(~10) 

-0.4840(-18) 

56. 

185.81993495 

3.0098488154 

-0. 12217 (-11) 

-0 . 8600 ( -19 ) 

57. 

71.17031124 

0.8742435461 

0.79800(-ll) 

0.2560(-18) 

58. 

317.03160513 

39.5236934204 

0 . 20929 ( -11) 

0. 3090( -18) 

59. 

212.45668639 

14.0505929610 

0.67903(-ll) 

0.2270(-18) 

60. 

343.34006144 

38.4321854270 

0 . 53277 (-11) 

0.4690( -18) 

61. 

73.24508049 

-  0.1113567166 

0.80957(-ll) 

0.3240(-18) 

62. 

303.05393222 

-26.1829392618 

-0 . 12252 ( -10) 

-0. 5440(-18) 

63. 

333.36324855 

-11.0407441456 

-0 . 80121 (-11) 

-0.3130(-18) 

64. 

81.17975482 

26.0770315311 

0 . 15371 ( -10) 

0.6360(-18) 

65. 

200.99943894 

36.6836983347 

-0. 10632(-10) 

-0 . 4300( -19) 

66. 

55.53089459 

49.7486910329 

-0.37261(-11) 

0.2520(-18) 

67. 

238.76514270 

12.9590849676 

0.10025(-10) 

0 . 3870( -18 ) 

68. 

339.19052294 

40.4033859525 

0 . 50962 (-11) 

0.3330(-18) 

69. 

50.32820349 

65.5477710863 

0. 19024 ( -10) 

0. 9910( -18) 

70. 

> 

0.0 

0.0 

0  0 

0.0 

71. 

0.0 

0.0 

J.O 

0.0 

72. 

0.0 

0.0 

0.0 

0.0 

73. 

0.0 

0.0 

0.0 

0.0 

74. 

0.0 

0.0 

0.0 

0.0 

!  75. 

1 

1 _ 

0.0 

0.0 

0.0 

0.0 

Note:  Arguments  numbered  70  through  106  are  all  zero  for  this  (B1950.0) 
version  of  the  transformation. 
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Table  2-3.  Coefficients  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 


a .  +  3  •  T 

A.  A, 

Y  .  +  6  .  T 

A,  A. 

-0.477766  (-02)  -  0.4838  (-05)  T 

0.255625  (-02)  +  0.247  (-06)  T 

0.572777  (-04)  +  0.5555  (-08)  T 

-0.24861  (-04)  +  0.1398  (-07)  T 

0.127777  (-05) 

-0.666666  (-06) 

0.305555  (-06) 

-0.833333  (-07) 

0.277777  (-07) 

-0.833333  (-07) 

-0.555555  (-07) 

0.277777  (-07) 

0.277777  (-07) 

-0.366305  (-03)  -  0.4444  (-07)  T 

0.159333  (-03)  -  0.8611  (-07)  T 

0.396111  (-04)  -  0.9444  (-07)  T 

0.150000  (-05)  -  0.2777  (-08)  T 

-0.143611  (-04)  +  0.3335  (-07)  T 

0.622222  (-05)  -  .1667  (-07)  T 

0.602777  (-05)  -  0.1388  (-07)  T 

-0.263888  (-05)  +  .8333  (-08)  T 

0.358333  (-05)  +  0.2777  (-05)  T 

-0.194444  (-05) 

0.133333  (-05) 

0.277778  (-07) 

-0.611111  (-06) 

0.472222  (-06)  -  0.2777  (-08)  T 

-0.416666  (-06) 

0.250000  (-06) 

-0.444444  (-06)  +  0.2778  (-08)  T 

0.194444  (-06) 

-0.333333  (-06) 

0.166666  (-06) 

-0.166666  (-06) 

0.833333  (-07) 

-0.138888  (-06) 

0.833333  (-07) 

0.111111  (-06) 

-0.555555  (-07) 
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-t 

+ 

—i 

y.  +  6.  T 

A.  A 

23. 

0.111111  (-06) 

-0.555555  (-07) 

24. 

-0.111111  (-06) 

25. 

0.277777  (-07) 

26. 

0.277777  (-07) 

27. 

-0.277777  (-07) 

28. 

0.277777  (-07) 

29. 

0.277777  (-07) 

30. 

-0.277777  (-07) 

31. 

-0.631668  (-04)  -  0.5555  (-08)  T 

0.271388  (-04)  -  0.1388  (-07)  T 

32. 

0.197777  (-04)  -  0.2777  (-08)  T 

-0.194444  (-05) 

33. 

0.107222  (-04)  -  0.1111  (-07)  T 

-0.555555  (-05) 

34. 

-0.836111K-05) 

0.358333  (-05)  -  0.2777(-08)T 

35. 

-0.438888  (-05) 

-0.277777  (-07) 

36. 

0.341666  (-05) 

-0.147222  (-05) 

37. 

0.175000  (-05) 

-0.555555  (-07) 

38. 

0.175000  (-05)  +  0.2777  (-08)  T 

-0.916666  (-06) 

39. 

-0.161111  (-05)  -  0.2777  (-08)  T 

0.888888  (-06) 

40. 

-0.163888  (-05) 

0.722222  (-06) 

41. 

-0.141666  (-05) 

0.750000  (-06) 

42. 

-0.105555  (-05) 

0.444444  (-06) 

43. 

0.805555  (-06) 

-0.277777  (-07) 

44. 

0.805555  (-06) 

-0.333333  (-06) 

45. 

-0.861111  (-06) 

0.361111  (-06) 
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able  2-3.  Coefficients  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 
(Con ' t) 


r 


i 

a  • 

|  A. 

46. 

0.122222  (-06) 

47. 

0.583333  (-05) 

!  48. 

0.444444  (-06) 

49. 

-0.361111  (-06) 

50. 

-0.277777  (-06) 

Y  .  +  6  -  T 
'a.  a. 


-0.277777  (-07) 
-0.277777  (-06) 
-0.222222  (-06) 
0.194444  (-06) 
0.1388888  (-06) 


51. 

52. 

53. 

54. 

55. 


-0.194444  (-06) 
0.194444  (-06) 
j-0. 194444  (-06) 
j-0. 222222  (-06) 
0.166666  (-06) 


-0.833333  (-07) 
0.833333  (-07) 
0.833333  (-07) 


56.  ! 

0. 

166666 

(-06) 

57.  i 

j 

-0. 

166666 

(-06) 

58. 

-0. 

194444 

(-06) 

59.  | 

i 

0. 

166666 

(-06) 

60.  | 
j 

j 

-0. 

138888 

(-06) 

61. 

1  0. 

138888 

(-06) 

62.  j 

j-0. 

,138888 

(-06) 

63. 

-0. 

.111111 

(-06) 

64.  ' 

j  0. 

.111111 

(-06) 

65. 

-0. 

1 

j 

.mm 

(-06) 

66 .  j 

j 

j-0. 

.833333 

(-07) 

6/. 

0. 

.833333 

(-07) 

f 

gp.  ^  j 

-0. 

.833333 

(-07) 

69.  1 

!  -f) . 
1 

.833333 

(-07) 

70.  1 

!-n. 

.  5 55555 

(-07) 

-0.833333  (-07) 
0.833333  (-07) 
0.833333  (-07) 
-0.833333  (-07) 
0.833333  (-07) 


0.833333  (-07) 


0.277777  (-07) 
0.277777  (-07) 
0.277777  (-07) 


Table  2-3.  Coefficients  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 
(Con't) 


l 

£ 

+ 

— i 

y.  +  6.T 

A.  A. 

B 

-0.833333  (-07) 

0.277777  (-07) 

B 

-0.833333  (-07) 

0.277777  (-07) 

BfStl 

0.555555  (-07) 

-0.277777  (-07) 

74. 

-0.555555  (-07) 

0.277777  (-07) 

75. 

0.555555  (-07) 

-0.277777  (-07) 

76. 

-0.555555  (-07) 

0.277777  (-07) 

77. 

0.555555  (-07) 

78. 

0.555555  (-07) 

-0.277777  (-07) 

79. 

0.277777  (-07) 

-0.277777  (-07) 

80. 

-0.277777  (-07) 

81. 

0.277777  (-07) 

-0.277777  (-07) 

82. 

-0.555555  (-07) 

0.277777  (-07) 

83. 

-0.277777  (-07) 

84. 

0.277777  (-07) 

-0.277777  (-07) 

85. 

-0.277777  (-07) 

0.277777  (-07) 

86. 

-0.27/777  (-07) 

0.277777  (-07) 

87. 

0.277777  (-07) 

88. 

0.277777  (-07) 

89. 

0.277777  (-07) 

-0.277777  (-07) 

90. 

-0.277777  (-07) 

91. 

-0.277777  (-07) 

92. 

0.277777  (-07) 

93. 

0.277777  (-07) 

94. 

-0.277777  (-07) 

95. 

0.277777  (-07) 
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Table  2-3.  Coefficients  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 
(Con 1 1) 


i 

a  ■  +  B  •  T 

A.  A. 

—  ■■■  . -  -■ 

Y  .  +  6  -  T 

A.  A. 

96. 

0.277777  (-07) 

97. 

-0.277777  (-07) 

98. 

-0.277777  (-07) 

99. 

-0.277777  (-07) 

100. 

-0.277777  (-07) 

101. 

-0.277777  (-07) 

102. 

-0.277777  (-07) 

103. 

-0.277777  (-07) 

104. 

0.277777  (-07) 

105. 

-0.277777  (-07) 

..... 

106. 

0.277777  (-07) 

T  =  (JED  -  2451545.0)  /  36525 

ms-~* 

The  quantities  in  the  parentheses 

! 

j 

are  the  exponents  of  10. 

J 

j 

i 

1 

j 

i 

i 

I 

1 

l 

i _ _ _ 

.  •  ■ 
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Table  2-4.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 


i 

al 

bi 

cl 

d. 

-L 

1. 

125.04452220 

-0.0529537648 

0 . 15522 (-11) 

0 . 4560 ( -19) 

250.08904440 

-0.1059075296 

0. 31044( -11 ) 

0. 9120( -19) 

41.66237996 

0.2757608151 

-0. 17006 ( - 10) 

-0 . 5586 ( - 18 ) 

83.38214224 

-0.3287145800 

0 . 18559 ( - 10 ) 

0.6042(-18) 

166.70690216 

0.2228070503 

-0.1 5454 ( -10 ) 

-0.5130(-18) 

-160.41505000 

-0.1113564531 

0 . 80742 (-11) 

0.3249(-18) 

-279.16783004 

0.05304790896 

-0.85852(-12) 

0 . 9120( - 19 ) 

8. 

208.42666444 

-0.3816683448 

0.20111(-10) 

-0.6498( -18) 

9. 

-159.06786124 

1.9712947103 

0. 45340( -12) 

0.0 

10. 

357.527723330 

0.9856002831 

-0.12014(-12) 

-0.6840(-19) 

11. 

198.45986206 

2.9568949934 

0 . 33326( -12 ) 

-0. 6840( -19) 

12. 

-156.59558454 

0.9856944272 

0 . 57354 ( -12 ) 

0. 6840( -19) 

13. 

-284.11238344 

2.0242484751 

-0.1 0988 (-10) 

-0.4560(-19) 

14. 

-325.77476340 

1.7484876600 

0 . 15908 ( -10 ) 

0. 5130( -18) 

15. 

-63.54522384 

26.5646080096 

-0.86250( -11) 

-0 . 3420 ( -19) 

16. 

355.05544660 

1.9712005662 

-0.24028(-12) 

-0. 1368( -18) 

17. 

122.57224550 

0.9326465182 

0.14320(-11) 

-0.2288( -19) 

18. 

195.98758536 

3.9424952765 

0 . 21312 ( -12) 

-0. 1368(-18) 

19. 

-232.48320110 

-1.0385540479 

0 . 16723 (-11) 

0.1 140 ( - 18 ) 

20. 

90.81928560 

-1.8014414248 

-0.14355(-10) 

-0 . 4674 ( -18 ) 

21. 

-281.64010674 

1.0386481920 

-0.97866(-12) 

0 . 2280 (-19) 

22. 

-200.73024120 

1.6955338951 

0. 17460( -10) 

0 . 5586 ( -18 ) 

23. 

73.41533986 

3.0098487582 

-0.12189(-11) 

-0. 1 140 ( - 18 ) 

24. 

-162.88738170 

0.8742438300 

0 . 79540 (-11) 

0 . 2565 ( -18 ) 

25. 

31.75295990 

2.7340879431 

0 . 15788 ( - 10 ) 

0 . 4446 ( - 18 ) 

The  quantities  in  the  parentheses  are  the  exponents  of  10. 
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Table  2-4.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 
(Con't) 


t 

, 

a . 

A. 

b. 

A. 

c . 

A. 

d  ■ 

A. 

26. 

309.16440924 

10.9348369451 

0.10722 ( -10) 

0 . 5016 ( -18 ) 

27. 

-119.35472040 

22.4102976738 

-0.26293(-ll) 

0.3534(-18) 

28. 

247.61676770 

0.8796927534 

0 . 29842 (-11) 

0.2280(-19) 

29. 

287.93190390 

-0.9271975948 

-0 . 64018 ( -11) 

-0 . 2109 ( -18 ) 

30. 

-51.62918234 

3.0628025231 

-0 . 2771 1 ( -11 ) 

-0. 1596(-18) 

31. 

76.63286496 

26.3527929503 

-0.24162(-11) 

0.2 166 ( - 18 ) 

32. 

134.96298140 

13.0649929500 

0 . 65192 (-11) 

0. 3648( -18) 

33. 

311.58834276 

26.4057467151 

-0 . 39684 (-11) 

0. 1710( -19) 

34. 

211.59584636 

39.4177859003 

0.41030(-11) 

0 . 5814 ( -18 ) 

35. 

-100.73774480 

-11.3165052900 

0.93883(-ll) 

0 . 1482 ( -18) 

36. 

301.66988356 

13.2878000003 

-0.89354(-ll) 

-0 . 1482 ( -18 ) 

37. 

235.70072620 

24.3814982400 

-0 . 28696 (-11) 

0.2166(-18) 

38. 

260.00750360 

13.0120391851 

+0. 80714 ( -11 ) 

0.4104(-19) 

39. 

-  9.91839580 

-13.1179467148 

-0.49670( -11 ) 

-0 . 3192 ( -18 ) 

40. 

177.37060976 

37.6692982403 

-0 . 11805 ( - 10) 

0.6840(-19) 

41. 

86.55132416 

39.4707396651 

0.25508(-ll) 

0. 5358( -18) 

42. 

312.33359116 

50.7342911903 

-0 . 52885 ( -11) 

0.4332(-18) 

43. 

269.925962801 

26.1299859000 

+0 . 13038 ( -10) 

0 . 7296 ( -18 ) 

44. 

-24.10487984 

15.0362876603 

0. 69726( -11 ) 

0 . 3648( -18 ) 

45. 

346.55882776 

52.4827788503 

0. 10622 ( -10) 

0 . 9462 ( -18 ) 

46. 

186.54382056 

26.4587004800 

-0. 55206( -11 ) 

0 . 1254 ( -18 ) 

47. 

176.62536136 

13.3407537651 

-0.10487 ( -10) 

-0 . 1938 ( -18 ) 

48. 

225.78226700 

11.2635515251 

-0. 78366(-ll ) 

-0 . 1026 ( - 18 ) 

49. 

-335.69322260 

-11.3694590548 

0. 10941 ( -10) 

0.1 938 ( - 18 ) 

50. 

322.25205036 

63.8522379051 

-0.3 1872 ( -12 ) 

0 . 7  524 ( - 18 ) 
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Table  2-4.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 
(Con1 t) 


i 

a . 

A. 

b. 

a. 

c . 

A. 

d. 

A. 

51. 

-103.21002150 

-10.3309050069 

0. 92687 ( -11 ) 

0.7980(-19) 

52. 

74.16058826 

27.3383932334 

-0.25363(-ll) 

0 . 1482 ( - 18 ) 

53. 

79.10514166 

25.3671926672 

-0.22960(-ll) 

0.2580(-18) 

54. 

87.29657256 

63.7992841403 

0 . 12334 (-11) 

0 . 7  980 ( - 18 ) 

55. 

10.66370760 

37.4464911900 

0 . 36496 (-11) 

0 . 5814 ( -18 ) 

56. 

110.85810156 

28.1012806103 

0 . 13491 ( -10 ) 

0 . 7296 ( - 18 ) 

57. 

0.74524840 

24.3285444751 

-0.13174(-11) 

0.2622(-18) 

58. 

187.28906896 

50.7872449551 

-0.68380(-11) 

0 . 3876 ( -18 ) 

59. 

-149.14940204 

15.0892414251 

0.54204(-ll) 

0.3192(-18) 

60. 

-110.65620400 

-24.4344520048 

0. 44218( -11 ) 

-0.1710(-18) 

61. 

-222.56474190 

12.0793926669 

0 . 66394 (-11) 

0.4332(-18) 

62. 

221.51430556 

52.5357326151 

0 . 90701 ( -11 ) 

0 . 9006 ( - 18 ) 

63. 

-238.17300290 

-23.3458979569 

0 . 27494 (-11) 

-0.2850(-18) 

64. 

-51.58083916 

-13.3937075530 

0 . 12039 ( -11 ) 

0 . 2394 ( -18 ) 

65. 

297.85036310 

12.1907491200 

-0 . 14348 (-11) 

0. 1083(-18) 

66. 

132.49070470 

14.0505932331 

0 . 63991 (-11) 

0. 2964(-18) 

67. 

321.50680196 

39.5236934300 

0 . 99868 (-11) 

0.4902(-18) 

68. 

214.06812306 

38.4321856172 

0 . 42232 (-11) 

0.6498(-18) 

69. 

179.84288646 

36.6836979512 

-0.11 684 ( -10 ) 

0. 1368 ( - 18 ) 

70. 

-144.88144060 

-26.1829396648 

-0 . 1 1486 ( -11 ) 

-0 . 6340 ( - 18 ) 

71. 

121.52180916 

65.5477718003 

0. 17141 ( -10) 

0. 131 1 ( -17) 

72. 

314.80586786 

49.7486909072 

-0.51 656 ( -11) 

0. 5016( -18) 

73. 

209.12356966 

40.4033861834 

0 . 39829( -11) 

0. 5130( -18) 

74. 

-59.07536484 

-11.0407444748 

-0. 76180( -11 ) 

-0.4104(-18) 

75. 

34.97048500 

26.0770321351 

0. 14590(-10) 

0 . 7752 ( - 18 ) 
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Table  2-4.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 
(Con't) 


A. 

a . 

A. 

b. 

A. 

d. 

A. 

76. 

25.05202580 

12.9590854203 

0 . 96236 (-11) 

0.4560( -18) 

77. 

44.88894420 

39.1949788500 

+0. 19557 (-10) 

0 . 1094 (-17) 

78. 

14.48322806 

38.5435420703 

-0.38510(-11) 

0. 3249( -18) 

79. 

115.12606300 

-13.1709004796 

-0.34148(-11) 

-0 . 2736 ( -18 ) 

80. 

-336.4384710 

-35.6980035300 

0 . 12258 ( - 10 ) 

-0. 6840( -19) 

81. 

42.40762836 

24.6043052903 

-0.18324(-10) 

-0.2964(-18) 

82. 

53.07133596 

62.0507964803 

-0.14674(-10) 

0. 2850( -18) 

83. 

-201.47548960 

-22.6330105800 

0 . 18777 ( -10) 

0. 2964( -18) 

84. 

333.42284346 

16.0218879434 

0.68525(-ll) 

0.2964(-18) 

85. 

322.25205036 

63.8522379051 

-0  -  31872 ( -12 ) 

0.7524(-18) 

86. 

278.10835456 

48.9855374251 

-0.21193(-10) 

-0. 7980( -19) 

87. 

128.21370412 

39.7465004803 

-0 . 44456 (-10) 

-0.2280(-19) 

88. 

-98.26546810 

-12.3021055731 

0 . 95090 (-11) 

0.21 66 ( -18 ) 

89. 

-14.18642064 

28.1542343751 

0. 11939( -10) 

0.6840(-18) 

90. 

222.25955396 

76.8642770903 

0 . 77527 (-11) 

0 . 1 162 (-17) 

91. 

135.70822980 

37.3935374251 

0. 52018 (-11) 

0 . 6270 ( -18) 

92. 

27.47595932 

28.4299951903 

-0.50672(-ll) 

0 . 1254 ( -18 ) 

93. 

245.82108296 

41.1662735603 

0.20011(-10) 

0 . 1094 ( - 17 ) 

94. 

-274.19392424 

15.1421951900 

0. 38682 (-11) 

0 . 2736 ( -18 ) 

95. 

309.11606606 

27.3913469982 

-0 . 40885 (-11) 

0 . 1026 ( - 18 ) 

96. 

228.25454370 

10.2779512420 

-0 . 77 1 65 ( -11) 

-0. 3420( -19) 

97. 

-61.49929836 

-26.5116542448 

0 . 70728 (-11) 

-0. 7980( -19) 

98. 

138.78250186 

14.1620438303 

-0.98140(-12) 

0 . 1 083 ( - 18 ) 

99. 

233.22844950 

25.3670985231 

-0.29897(-ll) 

0. 1482(-18) 

100. 

-287.28156536 

-37.7752057700 

0. 14909 ( -10 ) 

0 . 2280 ( - 19 ) 
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Table  2-4.  Arguments  of  the  Nutation  Sine  and  Cosine  Series  (J2000.0) 
(Con't) 


i 

a  ■ 

b. 

- c 

c  . 

-c 

d. 

g 

-45.93938054 

25.4201464320 

-0 . 38482 (-11) 

0 . 2394 ( -18 ) 

Bill 

21.83450070 

-10.3838587717 

0.10820(-10) 

0 . 1254 ( -18 ) 

103. 

184.11988704 

10.9877907100 

0 . 91702 ( -11 ) 

0 . 4560 ( - 18 ) 

104. 

145.62668900 

50.5114841400 

0 . 10168 ( -10) 

0.9462(-18) 

105. 

188.03431736 

75.1157894303 

-0.81554( -11) 

0.6498(-18) 

106. 

295.37808640 

13.1763494031 

-0 . 15549 ( -11) 

0.3990(-19) 
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4  THE  TRUE  OF  DATE  INERTIAL  REFERENCE  SYSTEM  AND  THE  C  TRANSFORMATION  (con* t) 


*e  x"  is  the  vector  in  the  true  of  date  system.  Due  to  the  orthogonality 
one  may  write  for  the  inverse  transformation: 


(13) 


EQUATORIAL  EARTH-FIXED  COORDINATE  SYSTEMS 
erence  systems  (iv),  (v),  and  (vi)  above  are  called  equatorial  earth¬ 
ed  systems.  Their  origins  lie  within  the  earth  and  they  utilize  the 
restrial  poles  and  equator  as  fundamental  components.  These  coordinate 
terns  are  not  inertial,  since  they  are  fixed  and  rotate  with  the  earth, 
h  of  these  three  reference  systems  and  associated  transformation  is 
cribed  in  the  following  subsections. 

.1  THE  INSTANTANEOUS  EARTH-FIXED  REFERENCE  SYSTEM  AND  THE  B  TRANSFORMATION 
instantaneous  earth-fixed  system  is  an  earth  centered  coordinate  set 
ch  has  its  (e‘,  f  ,  g')  axes  defined  as  follows: 

e'  is  the  unit  vector  in  the  true  equatorial  plane  in  the  direction 
of  the  Greenwich  meridian; 

■V  is  the  unit  vector  along  the  earth's  true  axis  of  rotation  and 
is  positive  in  the  northern  hemisphere;  and 

f  is  the  unit  vector  which  completes  the  right-handed  orthogonal 
system  with  e'  and  g' . 

geometry  of  Figure  2-5  depicts  the  relationship  between  the  true  of 
e  inertial  system  and  the  instantaneous  earth-fixed  system. 
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3.2  SMTP  STATION  CALIBRATION  (con't) 


p/c  =  transmission  delay,  found  as  described  in  3.1, 


t ,  =  receiver  equipment  delay, 

d . 

J 

t  =  station  clock  time  at  receipt  of  NAVSAT  time  mark,  and 
rec 

tc  r  *,  =  clock  error  of  j th  station  clock. 

The  linear  fit  is  made  between  Atcfa  and  tsta^  values  of  k  for 

iLdjk  3k 

each  station  j  for  some  epoch  tn  .  If  the  station  clock  is  reset  during  the 

j 

fit  interval,  it  must  be  treated  the  same  way  as  described  for  a  NAVSAT  clock 
reset  for  purposes  of  obtaining  a  fit  and  computing  the  time  error  afterward. 

«sta.  ’  6tstani  +  +  bsta,  +  '«j- 

J  oj  J  J 


where 


station  clock  error  in  seconds  at  time  t. 


'staoj 


constant  term  from  linear  fit  for  station  clock  j, 


0  before  station  clock  reset 


reset  adjustment  after  clock  reset 
(cumulative  if  more  than  one  reset). 


first  power  term  from  linear  fit,  and 


=  epoch  for  linear  fit  for  station  j. 


The  above  calibration  permits  correction  of  the  station  clock  time  at  the 
end  of  each  doppler  observation  interval,  i.e. 


-45- 


TR  82-391 


NAV SAT 


so  . 

1 


-t 


b 


si 


t 


CLOCK  CALIBRATION  (con't)  j 

=  constant  term  from  fit  for  NAVSAT  .j  ,  ~  1 

V-j 

=  0  before  NAVSAT.  reset,  •  ;.--j 

=  reset  adjustment  after  reset  (cumulative  if  more  than 

one  reset) ,  j 

J 

4 

=  first  power  term  from  linear  fit,  and  .j 

=  epoch  for  1  i  near  fi  t .  •; 


SMJ  P  J  T  AT  I_0N__CALJ  3  RATION 

r  each  NAV SAT  clock  has  been  calibrated,  each  station  clock  is  calibrated 
g  a  similar  procedure.  The  difference  in  the  two  procedures  is  that  .it 
sed  to  correct  for  errors  in  NAVSAT  time. 


ation  clock  offset  is  round  for  each  NAVSAT  time  mark  for  each  station 


'Sujk  '  tsuok  '  tfiAVSATi  *  **« 


(3) 


'sta¬ 


lk 


station  clock  error  for  k^  observation  at  station  j. 


'sta 


station  clock  time  at  NAVSAT  transmission  time  for  k 


th 


jk  observation, 


'rec  ‘  p/c  *  ld.  tscj  ’ 


ti~e  of  rearest  multiple  of  120  seconds, 


1  K 


correction  computed  from,  (2)  for  NAVSAT.  at 

1  .VMVSAl 


1 
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3.1  NAVSAT  CLOCK  CALIBRATION  (con't) 

Atscik  =  tNAVSATi'tstak 

where 

At  =  NAVSAT  clock  error  for  k^h  observation  of  NAVSAT  i, 

SCi  k 

lstak  =  trec  '  p/c  "  td-  '  tscj  > 

•=  reference  station  clock  time  when  signal  k  transmitted, 
tNAVSATj  =  time  of  nearest  multiple  of  120  seconds, 

p/c  =  transmission  time, 

t^  =  receiver  equipment  delay  for  reference  station  receiver, 

t  =  reference  station  clock  time  of  receipt  of  time  mark  ,  and 

'  f-  h 

t,ri  =  clock  error  of  jr  station  clock. 

bCJ 

An  array  is  created  for  each  NAVSAT  satellite  which  contains  all  observations 
of  (1)  transmit  time  (NAVSAT  clock)  and  (2)  the  NAVSAT  time  error. 

If  a  satellite  (NAVSAT)  clock  is  reset  during  the  interval  of  collection  of 
data,  the  NAVSAT  errors  following  reset  are  shifted  by  the  negative  of  the 
reset  (/,tc)  to  conform  to  the  earlier  data  for  purposes  of  the  linear  fit. 

A  linear  fit  is  made  to  the  data  set  for  each  NAVSAT  satellite.  This  linear 
fit  provides  the  ability  to  compute  the  NAVSAT  clock  error  at  any  time  during 
the  fit  interval.  The  value  derived  from  the  computation  is 
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SECTION  3 

TIME  CALIBRATION  AND  CORRECTION 


ange  difference  data  comes  from  the  tracking  stations  tagged  with  station 
lock  time  at  the  end  of  the  doppler  count.  However,  because  of  errors  which 
ay  exist  in  station  clocks,  it  is  necessary  to  have  a  means  of  calibration 
f  each  station  clock  such  that  the  time  tags  may  be  adjusted  before  any 
valuation  of  the  range  difference  data  can  be  made. 

he  calibration  is  accomplished  through  use  of  accurate  reference  clocks 
n  a  subset  of  the  tracking  stations  and  the  NAVSAT  clocks.  The  NAVSAT  clocks 
(ill  first  be  calibrated  by  use  of  the  reference  station  clocks.  Then  the 
■tation  clocks  will  be  calibrated  from  the  adjusted  NAVSAT  clock  times.  A 
inear  fit  is  used  to  calibrate  each  NAVSAT  clock  such  that  fit  coefficients 
letermine  the  error  at  any  time  during  the  fit  span.  Similarly,  a  linear  fit 
is  used  to  calibrate  each  station  clock  for  all  useable  NAVSAT  time  marks, 
before  each  fit,  a  check  is  made  to  evaluate  the  self-consistency  of  the  data, 
’oints  not  within  reasonable  bounds  are  deleted  from  the  fit  process. 


5.1  NAVSAT  CLOCK  CALIBRATION 

iach  NAVSAT  clock  is  calibrated  approximately  once  per  day.  Since  NAVSAT 
:locks  emit  time  signals  at  two  minute  intervals,  the  calibration  of  the 
1AVSAT  clock  must  be  made  with  time  corrected  from  receipt  to  transmit.  This 
idjustment  is  accomplished  by  subtracting  he  propagation  delay  and  the  re- 
:eiver  equipment  delay  from  the  calibrated  station's  received  time.  The  propa- 
jation  delay  is  found  by  using  SPASUR  elements  to  obtain  the  NAVSAT  position 
ising  the  Brouwer-Lyddane  prediction  method.  Then  the  range  (p)  is  found  by 
:alculating  the  distance  from  the  satellite  to  the  calibrated  station.  The 
equipment  delay  is  a  function  of  the  station  and  is  kept  in  the  station  table. 


rhe  NAVSAT  clock  error  is  computed  for  each  NAVSAT  time  mark  for  each  reference 
ibservation  station  as  follows: 


-  transformation  of  the  station  position  from  the  earth-fixed  to  the  inertial  frame  for 
pcomputa t ion  is  performed  using  equations  (29)-(30)  of  Section  2  when  X  is  replaced  by 
t_tVe)  +  X},  where  w  is  the  earth  rotation  rate,  t  the  time  from  the  beginning  of  the 
,  and  tve  the  time  from  the  beginning  of  the  day  of  the  first  transit  of  the  vernal 
inox.  The  es,  fs,  and  gs  components  are  now  inertial  station  components. 
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FIGURE  2-8.  t  Reference  Frame  Geometry 
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2.3  TIME  OF  CLOSEST  APPROACH  (tCA)  COORDINATE  SYSTEM  (con't) 

where  pUca)  and  p(tCA)  are  the  position  and  velocity  of  the  satellite 
relative  to  the  observing  station  at  the  time  of  closest  approach.  The 
geometry  of  the  tCA  system  is  shown  in  Figure  2-8. 

A  similarity  transformation  is  used  in  the  point  and  cross  pass  editors  to 
transform  the  normal  equations  to  the  reference  frame  (see  Section  7). 
The  matrix  used  in  this  transformation  is: 


In  the  above  expression  x,  y,  and  z  are  unit  vectors  directed  along  the 
respective  axes  of  the  inertial  reference  system. 


(33) 
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2.2.4  THE  GEOGRAPHICAL  REFERENCE  SYSTEM 

Cartesian  coordinates  for  a  station  in  the  earth-fixed  reference  ellipsoid 
system  are: 


(30) 


and  ag  and  e  are  the  semi -major  axis  and  eccentricity  of  the  reference 
ellipsoid,  respectively.  The  station's  coordinates  in  the  instantaneous 
earth-fixed  or  any  of  the  inertial  systems  can  be  obtained  by  the  proper 
application  of  the  A,  B,  C,  and  D  transformations  as  discussed  in  the  pro¬ 
ceeding  sections. 

2.3  TIME  OF  CLOSEST  APPROACH  COORDINATE  SYSTEM 

The  t^  reference  system  is  used  extensively  in  the  point  and  cross  pass 
editors  (see  Section  7).  This  coordinate  system  is  formed  using  the  position 
and  velocity  vectors  of  the  satellite  relative  to  the  observing  station  at 
the  time  of  closest  approach,  the  computation  of  which  is  described  in  Section 
9.  The  axes  of  this  time  dependent  reference  frame  are  defined  as  follows: 


)  x  p(t 


(31) 
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FIGURE  2-7.  The  Geographical  Coordinate  System 
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2.2.3  UT1  -  UTC  AND  POLE  WANDER  PREDICTION  (con't) 


(28) 


where  t  is  the  time  in  days  and  A'.,  B'.,  C‘-  ,  and  D'-  ,  are  externally  supplied 

■L  -t  ■<-£ 

expansion  coefficients,  and  P^  are  periods  in  days.  The  AX^-  in  equation  (28) 
represents  the  pole  positions,  i.e.  Ap  when  i=  1  and  A q  when  i=2. 


2.2.4  THE  GEOGRAPHICAL  REFERENCE  SYSTEM 

The  geographical  reference  system  uses  two  angular  measures  (one  from  the 
Greenwich  meridian  and  one  from  the  equatorial  plane  of  the  earth-fixed  ref¬ 
erence  ellipsoid  system)  and  a  linear  measure  (of  height  above  tne  refer¬ 
ence  ellipsoid)  to  specify  the  location  of  an  observing  station  on  the  surface 
of  the  earth.  These  elements  of  location  are  defined  as: 

<p  =  geodetic  latitude  (the  angular  measure  between  the  equator 

and  the  normal  to  the  tangent  surface  at  the  station  location); 

X  =  longitude  (the  angular  measure  between  the  Greenwich  meridian 
and  the  station  meridian  as  measured  in  the  equatorial  plane); 
and 

h  =  geodetic  height  of  observing  station  above  the  reference 
ellipsoid. 

The  geometry  of  this  reference  system  is  depicted  in  Figure  2-7. 
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2.2.2  THE  EARTH-FIXED  REFERENCE  ELLIPSOID  COORDINATE  SYSTEM  AND  THE  A 
TRANSFORMATION  (con't) 

L  =  Ax  ■  *  A  B  X"  =  A  B  C  X’  =  ABCDX,  (24) 

where  Xg  is  the  vector  in  the  earth-fixed  reference  ellipsoid  system.  The 
inverse  transformation  may  be  written  as: 

V  =  ATXe  (25) 

or 

X  =  DT  CT  BT  AT  Xe  .  (26) 

2.2.3  UT1  -  UTC  AND  POLE  WANDER  PREDICTION 

PULSAR  uses  externally  computed  expansion  coefficients  to  generate  and 
predict  values  for  A p,  A q,  and  At  at  one  day  intervals.  Values  for  those 
quantities  at  intermediate  times  are  obtained  using  the  linear  interpola¬ 
tion  method  discussed  in  Section  9. 

The  UT1  -  UTC  predictions  are  made  using  the  expansion 

k 

at  =  AtBt+^[Cm  s1n(-v)  tDmcos  (t^)]  (27> 

m=l 

where  t  is  the  time  in  days  and  A,  B,  Cm>  and  Dm  are  externally  supplied 
expansion  coefficients. 

Pole  position  predictions  are  made  using  the  expansion 
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2.2.2  THE  EARTH-FIXED  REFERENCE  ELLIPSOID  COORDINATE  SYSTEM  AND  THE  A 
TRANSFORMATION 

For  reasons  not  completely  understood,  the  crust  of  the  earth  slips  slightly 
about  the  earth's  axis  of  rotation.  This  produces  an  effect  called  polar 

AAA 

wander,  which  is  the  displacementof  the  true  pole  of  the  (e‘,  f',  g ' )  system 
from  the  pole  of  the  earth-fixed  reference  ellipsoid  system.  This  is  a  geo¬ 
centric  coordinate  space  which  has  its  (e,  f,  g)  axes  defined  as  follows: 

A 

e  is  the  unit  vector  in  the  plane  of  the  equator  at  the  date  of 
observation  and  is  in  the  direction  of  the  Greenwich  meridian; 

A 

g  is  the  unit  vector  along  the  semiminor  axis  of  the  reference 
ellipsoid  and  is  positive  in  the  northern  hemisphere;  and 

A 

f  is  the  unit  vector  which  completes  the  right-handed  orthogonal 

/V  A 

system  with  e  and  g. 

The  geometry  of  polar  wander  is  shown  in  Figure  2-6.  The  displacements 
A p  and  A q  are  angular  quantities  measured  from  the  pole  in  the  direction 
of  Greenwich  and  perpendicular  to  this  direction  (positive  towards  the  west), 
respectively. 

The  A  matrix  corrects  for  polar  wander  by  transforming  from  the  instantan¬ 
eous  earth-fixed  system  to  the  reference  ellipsoid  system.  This  matrix  is 
given  by: 


An  arbitrary  vector  x  may  be  transformed  as  follows  using  the  A  matrix: 
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2.2.1  THE  INSTANTANEOUS  EARTH-FIXED  REFERENCE  SYSTEM  AND  THE  B  TRANSFORMATION 
( con ' t ) 

where  e  and  At \>  are  given  by  equations  (7),  (8),  (9)  and  (10).  H0  is  the 
mean  hour  angle  of  Greenwich  at  the  beginning  of  the  day  of  observation 
(0h  UT) 2 .  It  is  computed  from 

H0  =  mod(GMS,  360°), 

where  GMS  is  the  observed  Greenwich  mean  sidereal  angle: 


GMS ( 0 )  = 


100.000100269  +  0. 9856473459844(t-t  )  +  0.29015  x  10'  (t-t  )2, 
B1950.0  0  0 


100.460618333  +  0. 9856473662902 (t-tD)  +  0.29077  x  10"  ( t-t0)2  . 

J2000.0 


Using  the  B  matrix  one  may  write  for  an  arbitrary  vector  V: 


V  = 


BCD  V 


where  V  ,  is  the  vector  in  the  instantaneous  earth-fixed  reference  system. 
From  the  orthogonality  property  of  the  B  matrix  one  may  also  write: 


T 

$  V 


T  T  T  ■* 

D1  C1  B1  V  , 
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2.2.1  THE  INSTANTANEOUS  EARTH-FIXED  REFERENCE  SYSTEM  AND  THE  B  TRANSFORMATION 
(con't) 

Transformation  from  the  true  of  date  system  to  the  instantaneous  earth 
fixed  system  is  accomplished  through  application  of  the  B  matrix  given  by: 


cos  A  sin  A 

-sin  A  cos  A 


where  A  is  the  longitude  of  the  Greenwich  Meridian  from  the  true  vernal 
equinox  of  date  and  is  given  by: 

A  =  Hq  +  AH  +  u(tf  -  At), 

where  w  is  the  earth's  mean  sidereal  rotation  rate  determined  by  the  time 
lapse  between  successive  transits  of  the  mean  equinox  and  tf  is  time  of  obser¬ 
vation  from  the  beginning  of  the  day. 

At  is  the  UTC  -  UT1  correction  2, 3, 4, 5, 6  for  -j rregul ar  rotation  of  the 
earth.  Values  for  this  are  published  by  the  Bureau  International  de  l'Heure. 
The  AH  in  the  last  expression  is  the  equation  of  the  equinoxes  shown  on 
Figure  2-4  and  is  defined  as: 


AH  =  tan-1  £tan  A^  cosej  , 


Universal  Time  (UT)  is  the  local  mean  solar  time  of  the  Greenwich  meridian. 

UT0  is  Universal  Time  as  determined  from  observations. 

UT1  is  Universal  Time  obtained  by  correcting  UT0  for  polar  motion. 

UT2  is  Universal  Time  obtained  by  correcting  UT1  for  seasonal  variations 
in  the  Earth's  rotational  speed. 

UTC  is  coordinated  Universal  Time  which  is  an  approximation  to  UT2  and  is 
based  upon  the  atomic  standard.  This  time  is  broadcast  by  government  agencies. 
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3.2  SMTP  STATION  CALIBRATION  (con't) 

T  .  =  t  .  -At  . 

obs  sta .  sta.* 

J  j 

where 

t  =  time  on  input  observation  file  for  observation  from  station  j 
j  =  station  clock  time  at  end  of  doppler  observation  interval. 

The  station  clock  calibration  is  done  on  all  SMTP  stations,  including  the 
reference  stations.  All  geodetic  satellite  doppler  data  from  all  stations 
are  adjusted  using  the  above  procedures. 
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SECTION  4 

OBSERVATION  CORRECTIONS 
4.1  TROPOSPHERIC  CORRECTION 

The  method  used  for  tropospheric  correction  in  PULSAR  is  a  version  of  the 
Hopfield  model  which  has  been  modified  for  the  pt.-pose  of  simplifying  nu¬ 
merical  computation.  The  associated  integral  is  evaluated  in  closed  form, 
thus  precluding  integration  of  each  data  point.  The  basic  modification  made 
is  in  the  expression  for  the  refraction  term  N.  The  expression  used  is: 

N  =  n  -  1  =  Nq 

where 

n  =  index  of  refraction, 

r-j-  =  upper  bound  of  integration, 

rQ  =  station  distance  from  center  of  earth,  and 

r  =  variable  distance  from  center  of  earth  along  line  of  sight 
from  the  station  to  the  satellite. 

The  above  expression  applies  both  to  the  dry  term  and  the  wet  term.  The  N0 
and  rT  will  be  different  for  the  two  cases. 

The  integral  to  be  evaluated  is 

y*  *  dp 

(n-l)  ~  — _  (earth-fixed  parameters) ,  (2) 

rsta  (L  *  r  -  k2)h 
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4.1  TROPOSPHERIC  CORRECTION  (con't) 


E/t2 

For  the  wet  term,  N0w  =  (.373  *  10‘2)  K  and  rT  ,  the  effective  upper 

bound, is  given  by: 


rT  =  r  +  12.0  (kilometers) 
'  w  0 


In  the  above  expressions 


r  _  n  Q(-37 .2465  +  .213166T  -  .000256908T2) , 

E  =  H  e'  k  k 


H  =  relative  humidity,  percent  (default  =  50), 


P  =  pressure,  millibars  (default  =  980), 


Tk  =  T  +  273  ,  and 


T  =  temperature,  °celcius  (default  =  15) 


The  total  correction  to  range  is  the  sum  of  the  two  terms 


AfR  =  ARd  +  ARW  • 


The  actual  range  is  always  less  than  that  computed  from  signal  transit  time; 
so,  Afp  represents  the  amount  by  which  the  measured  range  is  excessive.  For 
range  difference  data,  this  correction  term  is  computed  at  both  ends  of  the 
observation  period  and  differenced  to  find  the  correction  to  range  difference. 
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4 . 2  ANTENNA  OFFSET  CORRECTION 

Since  satellite  transmitting  antennae  locations  are  normally  displaced 
from  the  center  of  mass,  quaternions  relating  body  axes  to  inertial  coor¬ 
dinates  can  be  used  to  make  the  corrections  that  result  from  such  displace¬ 
ments  to  the  computed  range  values. 


For  purposes  of  these  computations,  it  has  been  found  adequate  to  interpolate 
over  a  set  of  eight  four  element  quaternion  arrays  to  find  the  elements  of 
the  array  at  any  desired  time.  However,  the  interpolated  array  must  be  nor¬ 
malized  such  that  it  has  an  absolute  value  of  unity.  The  normalization  may 
be  done  by  dividing  each  element  by  the  square  root  of  the  sum  of  the  squares 
of  the  four  elements. 

A  quaternion  set  may  be  converted  to  a  direction  cosine  matrix.  The  elements 
of  the  direction  cosine  matrix  are 


q  ii  = 

2 

q22  ■ 

■  Q: 

12  + 

qi2  = 

2 

q2  + 

«3 

Q4) 

q13  = 

2 

Ml 

V 

Q4) 

q21  = 

2 

Mi 

Q2  - 

% 

q4) 

q22  = 

-Q 

f  ' 

h  q22 

-  ( 

h 2  ■ 

q23 

2 

m2 

^3 

«1 

Q4) 

q31 

2 

Q3  + 

Q2 

q4) 

q32 

2 

(Q2 

q3  - 

q4) 

q33 

-c 

>!*• 

-  V 

+  ( 

V 

(14) 


As  a  quaternion,  Q  has  the  structure  Q^i  +  Q2j  +  Qgk  +  Q^. 
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4.2  ANTENNA  OFFSET  CORRECTION  (con’t) 


The  Q  matrix  defines  the  orientation  of  a  set  of  orthonormal  vehicle  fixed 

a  A 

axes  [  b  ]  relative  to  an  inertial  reference  frame  [i  ]: 


1 

1“4 

-O 

1 - 

r  •  1 
7l 

b2 

t  Q  ] 

*2 

b3 

h 

(15) 


The  Q  matrix  is  used  to  convert  the  vector  originating  at  the  spacecraft 

center  of  mass  and  pointing  to  a  transmitting  antenna  to  a  vector  in  inertial 

coordinates,  Ar^.  Thus, 

i?A  =  l?,T  IV]bc  ’  (16) 


where  mbc  is  the  body  centered  vector  from  the  center  of  mass  of  the 
spacecraft  to  the  phase  center  of  the  antenna. 


be 


l 

L  -J 


(all  coordinates  body  centered) 


(17) 


cm 
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4.2  ANTENNA  OFFSET  CORRECTION  (con't) 


The  coordinates  of  the  antenna  phase  centers  will  be  constant  for 
a  particular  spacecraft. 


The  computed  vector 


r 

— 

p„ 

X 

X 

X 

Pw 

= 

Y 

Y 

y 

P2 

Z 

sat. 

Z 

-  _ 

CM 

—  — 

L  J  sta 


(18) 


must  be  adjusted  by  addition  of  Ar^  to  compute  the  range  from  the  station 
to  the  antenna  (see  Section  7). 


4.3  RELATIVISTIC  CORRECTIONS 


Range  measurements  and  range  difference  measurements  computed  classically 
from  signal  transit  times  are  in  error  because  of  relativistic  effects.  The 
expression  for  the  range  corrected  for  relativistic  effects  is1 


p  =  PU  -  apgr 


or 


tRt 


pu  + 


3p 


4 


'  —  »"(tRe)*\/(tRe),  (19) 


S-M 


1  Gibson,  L.  R.,  Periodic  Relativistic  Effects  in  Satellite  Tracking, 
NSWC/DL  TR  (In  preparation),  Dahlgren,  Virginia,  1981. 
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4.3  RELATIVISTIC  CORRECTIONS  (con't) 
where 

Pu  =  Uncorrected  range  computed  classically, 

t  =  station  clock  time  when  pulse  was  emitted  by  beacon, 
kg 

p  =  Earth's  gravitational  constant, 
rR  =  Station  distance  from  earth  center  of  mass, 
a j  =  Semi-major  axis  of  satellite  orbit, 

VR  =  Velocity  of  station  in  ECI  frame. 
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4.3  RELATIVISTIC  CORRECTIONS  (con't) 

The  first  term  is  the  uncorrected  range  difference,  the  second  term  varies 


constant  AtRg.  The  third  term  is  cyclic,  varying  with  the  eccentric  anomaly 
of  the  satellite.  It  would,  of  course,  be  zero  for  a  circular  orbit  and 


would  be  small  enough  to  be  negligible  for  many  near-circular  orbit  satel¬ 
lites.  However,  PULSAR  cannot  ignore  this  term.  So,  corrections  for  both 
the  second  and  third  terms  are  made. 

In  the  above  expression  the  uncorrected  range  difference  is 

Apu  =  JL  [-N  -  AfAtR  ]  ,  (21) 

where 

N  =  Number  of  doppler  counts  within  AtR, 

Af  =  fR  -  fT, 

fR  =  frequency  of  receiving  station, 

fy  =  frequency  of  transmitting  station,  and 

AtR  =  receiver  clock  interval  for  the  N  counts. 

Tests  have  been  made  to  determine  the  magnitude  of  the  correction  terms. 
Computations  were  made  using  satellite  orbital  elements  typical  of  those 
supported  by  PULSAR  for  range  difference  with  an  observation  time  of  twenty 
seconds.  The  AtRe  term  was  of  the  order  of  two  meters  range  difference 

-y  , 

correction.  The  last  term,  involving  r  •  v*  varied  in  range 
difference  correction  over  a  fourteen  minute  interval  (approximately  a 
station  pass)  from  a  Ap  correction  of  -.14  meters  to  .23  meters. 
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4.4  CORRECTION  FOR  ROTATION  OF  CIRCULARLY  POLARIZED  ANTENNA 

As  is  the  case  with  many  missions,  e.g.  ,  space  shuttle  and  LANDSAT,  periodic 
adjustment  of  vehicle  attitude  provides  the  possibility  that  rotation  about 
the  antenna  can  produce  erroneous  results  in  the  range  difference  computations. 
A  procedure  to  compensate  for  the  induced  error  has  been  devised  and  involves 
calculating  the  angular  rotation  of  the  satellite  about  the  antenna  boresight 
and  computing  the  apparent  doppler  shift  produced  by  that  rotation. 

Blanton2  has  shown  that  the  angular  rotation  can  be  computed  from  the 
orientation  quaternions  by 


„„  coAt 

COS  - 

A 

c 

= 

Qj(t)  Q2{t)  Q3(t)  Q4(t) 

Q^t+At) 

ojjSinuAt 
a)  2 

= 

Q4(t)  Q3(t)  -  Q2(t)  -  Qj(t) 

Q2(t+At) 

<*)2  sin/joAt 
co  2  ’ 

= 

-Q3(t)  Q4(t)  Q1(t)  -Q2(t) 

Q  3 ( t+A  t ) 

u)3sin<Mt 

co  2 

—  — 

= 

Q_(t)  -Q  (t)  q  (t)  -Q  (t) 

2  14  3 

Q4(t+At) 

=  [QSQ1  IQ] 

This  formulation  assumes  that  the  components  <o*  remain  constant  over  the 
interval  At.  The  accuracy  of  this  method  is  insensitive  to  the  size  of 
the  interval  At.  Thus,  it  is  only  necessary  to  use  the  quaternions  at  the 
beginning  and  end  of  an  observation  interval  so  that: 


2  Blanton,  J.  N.,  Memo  to  W.P.B.,  9/2/81 
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*•4  CORRECTION  FOR  ROTATION  OF  CIRCULARLY  POLARIZED  ANTENNA  (con't) 


w  =  angular  rotation  rate,  and 
ft  =  angle  of  rotation  in  radians. 


here 
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4.4  CORRECTION  FOR  ROTATION  OF  CIRCULARLY  POLARIZED  ANTENNA  (con ' t) 
Then 


ft/2  = 

cos-'*'  Pj  > 

(24) 

ft  = 

2  cos-1  ,  and 

(25) 

:  '■  : 

ft-j  - 

ft 

P,  ,  ft/si n( — )  . 

i+i  2 

(26) 

* 

The  three  components  of  the  rotation  make  it  possible  to  compute  the  rotation 
about  the  beacon  boresight,  ft^bs,  which  has  body  fixed  direction  cosines  a ^ 
a,,,  and  a^. 

Then  3 

ft,  .  =  \  ft-a. 

bbs  Z__,  1  1  (27) 

i  =  l 

Subsequent  references  to  ft  will  imply  fibbs- 

The  range  difference  error  for  a  single  frequency  beacon  is 

c  ft 

6AR  =  - • -  >  (28) 

f  2tt 

where  f  is  the  frequency  and  c  is  the  speed  of  light. 


The  actual  beacons,  however,  use  frequency  pairs  to  compensate  for  ionospheric 
effects  (first  order).  The  derivation  which  follows  assumes  a  400  mHz,  150  mHz 
pair.  The  results  are  the  same  whether  processed  by  a  TRANET  receiver  or  by 
a  GEOCEIVER.  So,  the  development  will  be  presented  only  once  using  notation 
which  could  be  appropriate  to  either  type  processing.  Let  the  associated  offset 
frequencies  be  given  by 
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CORRECTION  FOR  ROTATION  OF  CIRCULARLY  POLARIZED  ANTENNA  (con’t) 


-fl  .  dl 


_ +  _ 


rate  of  rotation  about  the  transmitting  antenna.  Combining  the 
et  frequencies  to  eliminate  the  ionospheric  correction  term  involving 


Af'  =  Afj  -  XAf2 


(i-A2)  p+(i-a)  , 

c  2n 


x  =  f2/f1 


ne  sets 


f1  =  fx( 1-X2 ) 


55 

64  1 


Af'  =  -_  p+  (1-X)  — 
c  2tt 


p  =  —  -Af'  +  (1-X)  JL 
f '  2tt 


=  f1'  pdt  =  ±  J ^  -Af '  dt  +  —  (1-X) 
ti-l 
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4.4  CORRECTION  FOR  ROTATION  OF  CIRCULARLY  POLARIZED  ANTENNA  (con't) 


Now 


Af 


1 


), 


where 


fR  =  frequency  of  the  signal  received, 


fr  =  receiver  reference  frequency, 


also. 


Af„  = 


transmitter  frequency. 


), 


•  • 

and  N.j  and  are  the  doppler  count  rates  for  the  upper  and  lower  frequen¬ 
cies  respectively. 


The  signs  above  are  correct  for  the  transmitter  frequency  being  greater 
than  the  receiver  local  oscillator  frequency.  If  fr^  >  fs-,  then  the  sign  of 
N*  is  reversed. 

i 


Substituting  the  above  into  equation  (31)  gives 


5. 2. 1.1  The  Canonical  Equations  of  Motion  (con't) 


conjugate  momentum  P^  is  a  constant,  i.e.  P^  =  a£  ,  and 


H  =  H  (p  ,  p  , 

1  2 


£-r 


p 


A+r 


p 


3n-fe 


V  V  q£-l’  q£+l*  q3n-fe;  ^  * 

Thus  there  are  2(3n  -  k)  -2  equations  to  be  solved  and  the  problem  has  in 
fact  been  reduced  in  complexity.  Consequently,  the  canonical  formulation 
is  particularly  well  suited  for  dealing  with  problems  in  which  one  or  more 
of  the  coordinates  is  ignorable.  It  is  clear  that  the  simplest  solution 
to  a  dynamical  problem  would  result  if  the  problem  could  be  formulated  in 
such  a  way  that  all  coordinates  were  ignorable. 


5.2. 1.2  Hamilton-Jacobi  Theory 

Because  of  the  simplification  of  the  Hamiltonian  by  ignorable  coordinates, 
it  is  desirable  to  seek  coordinate  transformations  which  would  make  all 
the  coordinates  ignorable  thereby  making  the  solution  trivial.  Thus  given 
equations  (27),  one  wishes  to  change  from  the  3n-fc  variables  (q-,p<)  to 
3n-fe  new  variables  (Q.y ,  Py )  in  such  a  way  that  it  will  be  easy  to  determine 
the  new  equations  of  motion  and  at  the  same  time  preserve  the  canonical  form 
of  equations  (27);  i.e.  one  seeks  a  canonical  transformation  which  produces 
a  new  Hamiltonian  H'  so  that 

3H'  .  3fT 

0.  =  -  ;  p.  -  -  -  (j  --  1 ,  1,  ...  3 n-k), 

J  3Py  S  3Qy 


where  h'  is  easily  related  to  H  ^nd  is  expressed  as  a  function  of  Q-  and  ?■ 

J  J 

According  to  a  theorem  devised  by  Carl  Gustav  Jacob  Jacobi,  the  desired 
transformations  between  (qy,  j.-y )  and  (0-,  ?■)  which  give  rise  to  the  canoni 
cal  equations  (29)  are 
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ROUWER  METHOD  (con't) 

s  in  a  more  analytically  tractable  problem  since  two  of  the  trans¬ 
coordinates  of  position  become  constants  of  the  motion. 

f  outline  of  the  underlying  theory  is  presented  in  section  5.2.1  and 
suits  of  the  Brouwer  approach  are  summarized  in  section  5.2.2.  It 
be  mentioned  that  the  method  used  in  PULSAR  is  that  of  a  Lyddane 
ed  Brouwer  orbit  generator  which  has  been  modified  to  cope  with  small 
1  eccentricities. 

UNDERLYING  THEORY 

onian  dynamics  and  the  application  of  perturbation  theory  are  the 
elements  used  in  the  development  of  the  Brouwer  theory.  In  the  fol- 
subsections  the  canonical  equations  of  motion,  Hamil ton-Jacobi  theory, 
e  Von  Zeipel  method  of  finding  generating  functions  will  be  briefly 
sed. 


1  The  Canonical  Equations  of  Motion 

ler  the  system  of  n  particles  subject  to  k  constraints  which  limit  the 
i's  freedom  of  motion.  The  Hamilton  equations  of  motion  which  describe 
stem  are  given  by 


•  3H  .  3H  ,  ,  .  . 

;.  =  — —  ;  P.  = - (j  =  1,  2,  ...,  3n  -  fe),  (27) 

J  J  3q  • 


H  =  H(p.;  q  t)  is  the  Hamiltonian  for  the  system  and  q-  and  p- 

j  j  J  J 

ie  associated  3n  -  k  generalized  coordinates  and  momenta,  respecti- 
Because  of  the  symmetrical  appearance  of  equations  ( 27) ,  they  are 
nown  as  the  canonical  equations  of  motion. 

>uld  be  pointed  out  that  in  the  canonical  formulation,  if  q^  is  an 
ble  coordinate,  i.e.  it  does  not  appear  in  the  Hamiltonian,  then  its 
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5. 1.1. 3  Luni-Solar  Gravitational  Accelerations  (con't) 


Replacing  each  of  r$  and  r^1  by  the  accelerations  exerted  by  the  sun  and 
moon,  one  finds  that 


sun  s 


r  +  r  r 

s  moon  s  moon 


Treating  the  sun  and  moon  as  point  masses,  one  may  write 


-u.  r  • 
x.  x.s 


r  -  r  . 

x. 


r.  =  U.  - , 

*  *■  13 

r  3 


where  i  =  sun,  moon  and  y^  is  the  gravitational  constant  for  the  y 
celestial  body.  Substituting  these  into  equation  (23)  gives  the  desired 
resul t 


a  +  a  =  r  =  -y 
sun  moon  Msun 


-y 

r  -  r 


r  -  r  1 3  Ir  1 3 

sun1  |rsun! 


r  -  r 


r  -  r_ 


5.2  BROUWER  METHOD 

The  drag-tree  Brouwer  method  is  an  analytic  technique  which  considers  only 
geopotential  acceleration  for  generating  the  orbit  of  an  artificial  earth 
satellite  and  makes  use  of  the  theory  of  generating  functions  in  its  develop 
ment.  This  approach  casts  the  Hamiltonian  for  the  system  into  a  form  which 


SATELLITE 


s  =  SATELLITE 
®  =  EARTH 
s  =  SUN 
(  =  MOON 


FIGURE  5-1.  Luni-Solar  Perturbation  Geometry 
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5. 1.1.2  Acceleration  Due  to  Atmospheric  Drag  (con't) 
where  k  is  the  geocentric  altitude  of  the  satellite  given  by: 

a 

h  =  r - - -  .  (20) 

(  1  +  . -2  -  sin2*)35 

'  1  -  e2  ' 

In  the  last  equation  *  is  the  geocentric  latitude  and  the  quantities 
{i  =  1,  2,  . . . ,  5)  in  equation  (19)  are  data  base  constants. 


A  few  comments  are  in  order  concerning  the  value  of  the  satellite  drag  co¬ 
efficient.  In  general  this  coefficient's  value  is  dependent  upon  the  physical 
properties  of  the  satellite's  surface,  as  well  as  the  satellite's  shape  and 
angle  of  attack  to  the  direction  of  motion,  so  that  satellites  which  operate 
with  continuously  varying  attitudes  should  exhibit  time  dependent  values  of 
Cq.  However  it  has  been  shown  that  CQ ~  2.2  for  satellites  with  orbital 
parameters  like  those  supported  by  PULSAR2.  Thus  PULSAR  models  Cp  (or  more 
precisely  CpA/W)  as  a  constant  over  a  trajectory  segment. 

5. 1.1.3  Luni-Solar  Gravitational  Accelerations 

To  obtain  expressions  for  the  gravitational  accelerations  exerted  by  the  sun 
and  moon  upon  the  satellite,  consider  the  geometry  depicted  by  Figure  5-1, 
where  0  is  some  fixed  origin  relative  to  the  earth.  Since 


then 

*  ■  V  -  V  •  <22> 


2  King-Hele,  D.,  Theory  of  Satellite  Orbits  in  an  Atmosphere,  pp.  12-17, 
Butterworth  &  Co.,  1964. 
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5. 1.1. 2  Acceleration  Due  to  Atmospheric  Drag  (con't) 

W  =  mg  =  weight  of  satellite; 

m  =  mass  of  satellite; 

g  =  gravitational  acceleration  at  the  surface  of  the  earth; 
v  =  |v|  =  speed  of  the  satellite  relative  to  the  atmosphere;  and 
p  =  local  atmospheric  density  at  the  satellite's  location. 

The  velocity  of  the  satellite  relative  to  the  rotating  atmosphere  is: 

v  =  r  -  u  x  r  ,  (17) 


where  r  and  r  are  the  satellite's  inertial  position  and  velocity,  respec¬ 
tively.  The  earth's  angular  velocity  in  the  basic  inertial  reference  frame, 
Z,  is  given  by: 


1 

(° ' 

1  o  ^ 

1  o  ' 

w  =  (BCD)"1 

1 

0 

=  (CD)_1B-1 

0 

1 

=  (CD)"1 

0 

( 5 1 

r  / 

l  V 

where  <3  is  the  earth's  mean  sidereal  rotation  rate  and  BCD  are  the  trans¬ 
formation  matrices  of  SECTION  2. 


The  atmospheric  density  used  in  equation  (16)  is  computed  using  a  simple 
density  model  which  considers  only  the  density  variation  with  altitude  above 
the  earth's  surface.  Specifically  the  atmospheric  density  is  computed  from 


(Vh-c  -  (c  liz+  c  (i  •  c  i1] 

e  L  1  2  3  4  5  j 


(19) 
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5. 1.1.1  Non-Central  Geopotential  Acceleration  (con't) 

To  compute  aE  PULSAR  systematically  evaluates  the  and  functions  by 

using  the  (n,  m)  s  (0,  0),  (1,  0)  values  of  these  functions  as  starting 

values  and  computing  (j^and  Vp  f°r  n  =  2,  3,  ....  nmax  from  the  horizontal 

stepping  recurrence  relations.  The  diagonal  stepping  relations  are  then 

applied  to  compute  (J*>  V*  and  the  (j\  V*  f°r  n  =  2,  3,  n  are  calcu- 

11  n  n  max 

lated  using  again  the  horizontal  stepping  relations.  This  procedure  is  repeated 

until  m  =  m  <n  .  The  acceleration  ac  is  then  obtained  by  application  of 
max  -  max  E 

equations  (8),  (13),  (14),  and  (15). 

5. 1.1. 2  Acceleration  Due  to  Atmospheric  Drag 

The  acceleration  arising  from  atmospheric  drag  as  used  by  PULSAR'S  force  model 
is  assumed  to  be  that  of  neutral  particle  drag  opposite  the  direction  of  motion 
in  an  earth  fixed,  wind-free  atmosphere.  Aerodynamic  lift  forces  are  ignored. 
Drag  acceleration  is  very  difficult  to  determine  accurately,  due  to  gross 
uncertainties  in  the  temporal  behavior  of  the  earth's  atmosphere  and  the 
value  of  the  satellite's  drag  coefficient.  Additional  inaccuracy  is  intro¬ 
duced  if  the  satellite  is  nonspherical  and  its  frontal  area  is  continually 
changing. 

The  drag  acceleration  acting  on  a  satellite  due  to  the  earth's  atmosphere 
as  modeled  in  PULSAR  is  given  by: 


1 

I 

a 

a 

i 


>  1 


an  =  -  ^Dgpvv 


where 


cda 

D  =  u  /W 


(16) 


Cp  =  the  dimensionless  drag  coefficient  for  the  satellite; 

A  =  cross-sectional  area  of  the  satellite  normal  to  the  direction 
of  motion; 
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5. 1.1.1  Non-Central  Geopotential  Acceleration  (con't) 

Since  the  integration  of  the  satellite  equations  of  motion  is  performed 
in  the  basic  inertial  reference  frame,  the  earth-fixed  acceleration  given 
by  equation  (6)  must  be  transformed  to  the  inertial  acceleration  by  applying 
the  BCD  transformations  of  SECTION  2  to  equation  (6).  Thus 


ar  =  (BCD) 


- 1 


vev 


=  (BCD) 


ES 

n=0  m=0 
(n/1) 


c  vQ 

nm  e 


tfj1  +  S 


m 


V  V 
nm  e  n 


(8) 


It  should  be  mentioned  that  summing  the  index  n  to  infinity  is  not  compu¬ 
tationally  realizible.  In  practice  the  summation  must  be  truncated  at 
n  =  n  ,  where  nm3v  is  a  user  selectable  integer  less  than  or  equal  to  30 

(of  course,  ar  is  more  accurately  represented  when  large  values  of  n 

t  ITlaX 

are  used) . 


For  the  sake  of  computational  convenience,  the  well  known  recurrence 

relations  for  the  associated  Legendre  functions  and  their  derivatives1  have 

been  used  to  obtain  recurrence  relations  for  l/",  Vm,  V  Um,  and  v  Vm  .  In 

n  n  e  n  e  n 

particular,  it  is  found  that: 


u 


m 

n+1 


,  ,m 

VI 


r(n-m+l) 


(2n+l)cos0  Um-  (n+m) 
n 


r(n-m+l) 


(2n+l)cos0 


(9) 

horizontal 

stepping, 

(10) 


1  Arfken,  G.,  Mathematical  Methods  for  Physicists,  p.  560,  Academic  Press, 
1970. 
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5. 1.1.1  Non-Central  Geopotential  Acceleration  (con't) 


where 


oo  11 

n=0  m=0 
(n/1) 


+  S_ 


nm 


Z!e 

-n+1 


P™  (cose)  cos  m\ 


n 

Pm  (cose)  sin  mA 
rn+l  n 


(3) 


(4) 


(5) 


The  acceleration  in  the  earth-fixed  reference  frame  is  obtained  by  taking 

the  corresponding  gradient  of  the  potential  in  equation  (3).  Since  Cnm 

and  S  are  constants,  one  obtains 
nm 


vV  =  V  V  [c  V  Um  +  S  V  Vml  ,  (6) 

e  /  v  ;  I  nm  e  n  nm  e  n  j 

n=0  m=0 
(ni*l) 

*>■ 

where  is  the  gradient  operator  in  the  earth-fixed  frame  and  is  given  by 


3 


In  this  expression,  the  e.  are  unit  vectors  along  the'  orthogonal  Cartesian 

■-C 

x-  axes  of  the  earth-fixed  reference  frame. 

t 
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5.1.1  FORCE  MODEL  (con't) 

amoon  =  acceleration  due  to  the  lunar  gravitational  field;  and 

aS(Jn  =  acceleration  due  to  the  solar  gravitational  field. 

The  following  sections  describe  the  manner  in  which  each  of  the  terms  on 

the  right  hand  side  of  equation  (l)  is  evaluated  in  PULSAR. 

5. 1.1.1  Non-Central  Geopotential  Acceleration 

The  earth's  gravitational  potential  in  an  earth-fixed  reference  frame  is 
given  by  the  spherical  harmonic  expression 

V"  ■  “  £  h  fe  p” (cos0)  lc™ cos  mA  +  s™ sin  “ )  1  •  (2) 

v  n=0  m^O  lr"  1  J 

where 

p  =  earth's  gravitational  constant; 
ae  =  earth's  semi -major  axis; 

r  =  magnitude  of  the  geocentric  radius  vector  to  the  satellite; 

P™(cose)  =  associated  Legendre  function; 

0  =  colatitude  of  sub-satellite  point; 

X  =  east  longitude  of  the  sub-satellite  point;  and 

Cnm,  Snm  =  gravitational  harmonic  expansion  coefficients. 

It  should  be  noted  that  with  the  geopotential  expressed  in  this  form 
CQ0  =  1;  and  C^q  =  =  0,  i.e.  n  f  1.  For  computational  con¬ 

venience,  equation  (2)  may  be  recast  into  the  form 
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5.1  COWELL  METHOD  (con't) 

significant  perturbations)  and  their  step-by-step  numerical  integration 
to  obtain  r(t)  and  r(t).  In  general  the  method  consists  of  the  following 
basic  steps: 

(i)  given  r  and  r  at  some  time  t^; 

( i i )  calculate  the  total  acceleration  acting  upon  the  satellite 
at  t^; 

(iii)  numerically  integrate  to  get  r  at  the  next  time  point  ^+1 ; 

(iv)  numerically  integrate  again  to  get  the  position  at  t^+j; 

(v)  with  the  new  r  and  r,  repeat  the  above  steps  for  the  next 
time  point,  t^+2- 


The  procedure  for  initializing  the  numerical  integration  and  the  compu¬ 
tational  techniques  employed  in  the  numerical  integration  process  are  devel¬ 
oped  in  Section  9.  A  discussion  of  the  satellite  equations  of  motion  as 
applied  to  the  orbital  force  environment  of  satellites  serviced  by  the  PULSAR 
system  is  given  below. 


5.1.1  FORCE  MODEL 

The  force  environment  for  the  satellites  serviced  by  the  PULSAR  system  can 
be  adequately  represented  by  the  acceleration  function 


r(r,f,t) 


=  a 


E  +  aD  +  amoon  +  asun 


(1) 


where 


a£  =  acceleration  due  to  the  gravitational  field  of  the  earth; 
ap  =  acceleration  due  to  atmospheric  drag; 
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SECTION  5 

EPHEMERIS  GENERATION 

In  the  context  of  the  PULSAR  data  editing  system,  ephemeris  generation  is 
the  process  of  calculating  positions  and  velocities  for  earth  satellites 
at  specified  times.  Although  this  ephemeris  data  need  not  be  of  superior 
geodetic  quality,  it  provides  for  the  editing  process  with  sufficient 
accuracy : 

(i)  computed  satellite  observations  used  for  state  estimation 
and  dynamic  parameter  improvement  (section  7); 

(ii)  time  histories  of  satellite  positions  and  velocities;  and 

(iii)  predicted  values  of  future  satellite  positions  and  velocities. 

Two  techniques  are  used  to  generate  earth  satellite  ephemerides  in  the  PULSAR 
system,  the  well  known  Cowell  and  Brouwer  methods.  The  Cowell  method 
is  a  special  perturbation  theory  requiring  the  numerical  integration  of  the 
satellite's  equations  of  motion  to  obtain  time-dependent  values  of  position 
and  velocity.  It  is  used  to  generate  the  satellite  trajectory  which  initiates 
the  point  editing  process.  The  Brouwer  method  is  a  general  perturbation 
theory  which  provides  the  earth  satellite's  ephemeris  through  an  analytic 
integration  of  a  modified  (and  analytically  tractable)  formulation  of  the 
satellite  equations  of  motion.  This  technique  is  used  in  time  correction  com¬ 
putations  of  the  time  calibration  procedure.  Both  of  these  basic  methods  are 
discussed  in  the  following  subsections. 

5.1  COWELL  METHOD 

The  Cowell  method  is  the  most  straightforward  for  determining  the  position 
r(t)  and  velocity  r(t)  for  the  perturbed  satellite  trajectory.  Its  appli¬ 
cation  involves  the  development  of  the  satellite's  equations  of  motion 
which  is  consistent  with  its  force  environment  (taking  care  to  include  all 
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4.4  CORRECTION  FOR  ROTATION  OF  CIRCULARLY  POLARIZED  ANTENNA  (con* t) 
Rearranging  and  integrating  yields 


AR- 


l 


^1  ”  ^2  c 

Tr774vvv(t; 


*1-1 }  + 


c 

(1  +A)f1 


n 

2tt 


The  contribution  due  to  rotation  is  the  last  term,  i.e. 


6AR- 

i 


8 

11 


C  Q 


However,  the  expression  for  ft  assumes  right  hand  circular  polarization 
With  left  hand  circular  polarization  the  sign  is  reversed  giving 
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5. 2. 1.2  Hamil ton-Jacobi  Theory  (con't) 


3S(Pj;  qj;  t) 


3 SlPjj  qj ;  T) 


(j  =  1,  2,  ...3n-fe),  (30) 


where  S ( py ;  qj ;  t)  is  called  the  generating  function.  The  new  Hamiltonian 
is  related  to  the  old  throuah  the  transformation 


H’  =  H 


5. 2. 1.3  Von  Zei pel's  Method 

In  this  and  the  remaining  subsections  of  section  5,  the  analytical  develop¬ 
ment  will  be  performed  in  terms  of  the  six  Delaunay  variables  defined  within 
the  context  of  PULSAR  by: 


-  yrr 


G  =  L  Vi  -  e: 


=  G  cosT. 


l  =  mean  anomaly 


g  =  argument  of  perigee 


h  -  right  ascension  of  the  ascending  node, 


where  y  is  the  earth's  gravitational  constant,  and  a,  e,  i  are  the  orbital 
semi-major  axis,  eccentricity  and  inclination,  respectively.  The  Delaunay 
variables  provide  the  following  system  of  canonical  equations  of  motion: 
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5. 2. 1.3  Von  Zeipel 1 s  Method  (con't) 

where  H  is  a  function  of  the  six  Delaunay  variables  (note  that  in  this 
context  the  variables  L,  G,  H  are  the  three  variables  of  momenta  and 
l ,  g,  h  are  the  three  variables  of  position). 

The  Von  Zeipel  method  involves  solving  equations  (33)  by  obtaining  a  series 
of  variable  transformations  with  the  aid  of  a  generating  function.  In  parti¬ 
cular,  it  is  seen  from  the  previous  section  that  if  one  considers  a  new  set 
of  variables  L‘,  G',  H',  l',  g' ,  k1  and  a  generating  function  S(L',  G',  H', 
l  ,  g  ,  k  )  such  that 


L  = 


as 

Ji 


as 

9L  1 


as 

3G 1 


as  as 

h’  =  -  ;  and  -  =  0; 

3H'  at 


(34) 


then  the  new  system  of  equations  is  also  canonical  and 

H(L,  G,  H,  l,  g,  h )  =  H'{  L  * ,  G\  H',  V,  g',  h' ) .  (35) 

The  transformation  is  selected  such  that  H‘  is  independent  of  one  of  the 
position  variables  on  which  H  depends.  If  this  operation  is  repeated  sev¬ 
eral  times,  all  position  variables  will  be  eliminated  from  H  one  by  one,  and 
the  final  Hamiltonian  H"  will  be  a  function  only  of  L",  G",  and  H".  Thus 
in  the  equations 

a  h"  .  a  w"  .  a  h" 

l"  =  -  ;  G"  =  -  ;  H"  =  -  ;  (36  ) 

3 l"  dg"  3 h" 

the  right  hand  sides  will  vanish  and  L",  G",  H"  will  be  constants.  Placing 

those  values  into  the  right  hand  sides  of  the  canonical  equations  for  the 

•  •  • 

position  coordinates,  one  finds  that  £",g",  and  h"  are  constants,  so  that 
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£  ,  g  ,  h.  are  linear  functions  of  time.  Relations  of  the  type  given  by 
equations  (34)  are  then  used  to  return  to  expressions  for  the  original 
(unprimed)  Delaunay  variables. 

In  order  to  solve  the  problem  in  this  fashion,  one  must  be  able  to  determi 
the  generating  function  S.  To  do  this  using  the  Von  Zeipel  method,  one 
expands  the  generating  function  and  Hamiltonians  in  powers  of  some  small 
parameter  e.  Then 


H  -  H0  +  Hi  +.. .  ;  S  =  S0  +  Si  +  S2 


H'  =  H0'  +  H, '  +  H2'  +  ...; 


and  if  one  selects 


I'l  +  G'g  +  H 'h  , 


the  application  of  the  second  of  equations  (37)  to  equations  (34)  gives 
L 


L'  ♦  !!i  ♦  !!•  + 


U 


U 


and  similarly  for  G  and  H;  and 


V  =  l  + 


as, 

3L 1 


3L  ’ 


and  similarly  for  g  and  h.  If  one  performs  a  Taylor  expansion  of 

u  (  as  as  as  ,  ,  \  .  . 

H  I -  ,  —  ,  - ;  g,  k)  about 

'  H  3g  3  h  I 


3So  3So  ^  as 

— 2-  ,  and  — 2_  and 


3  l 


3  h 


equates  terms  of  equal  orders,it  is  found  to  second  order  that: 
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5.2. 1.3  Von  Zei pel's  Method  (con't) 


Ho  H0'  (zeroth  order) 


3 


+  H2  =  H2  (second  order). 


where  Hx  and  Hz  are  functions  of  the  primed  coordinates  and  (a2 ,  a2,  a3)  = 

(L‘,  G1,  H’)  and  (gls  B2,B3)  =  (£,  g,  h) .  It  should  be  noted  that  this 

expansion  can  be  continued  indefinitely  to  higher  orders  of  accuracy.  These 
equations  are  then  solved  for  Si,  S2,  etc.,  depending  upon  the  accuracy 
required,  and  the  desired  solutions  follow  from  equations  (39)  -  (40).  This 
procedure  is  repeated  to  eliminate  the  remaining  position  coordinates. 

5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY 

To  illustrate  the  theory  outlined  in  section  5.2.1,  consider  its  application 
to  the  main  problem  of  the  artificial  satellite  theory,  i.e.  the  geopotential 
is  assumed  to  be  of  the  form 

"V"  =  T  +  “X-  (1  •  3  si"2  •  (42) 

1  2  1  2 

where  fe2  =  —  J2  a  =  —  C  a.  For  the  sake  of  brevity,  only  the  first 
'  2  e  2  20  e 

order  theory  through  the  el imination  of  l  will  be  discussed. 

The  Hamiltonian  for  the  system  is  given  by 
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5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 


p  life 

H  =  -  +  — 2-  (1-3  sin2?). 

2  a  r3 


Using  spherical  trigonometry  and  equations  (32),  one  may  rewrite  the  last 
equation  as 


(43) 


H  =  H0  +  Hj  , 


where 


2L" 


(44) 


(45) 


and 


\i  kz 

H,  =  - — 

1  L6 


1  +  JL  JlL\  g3  +  /  3  .  3  H2 

2  2  G2 


2  2  G2  /  r3 


(46) 


d 

-  cos(2g  +  2f) 

r3 


In  the  last  expression  f  is  the  true  anomaly.  It  may  also  be  shown  that 


*1  -  lL  +  V  2  P ' 

r3  G3  J 

y=i 


L3 


f  C0S  J  £  =  q3  ’  U1 


+  5, 


(47) 


and 


g 


cos  (2g  +  2f)  = 


E 

i  =  -• 
j  t  o 


cos  (2g  +  j£)  = 


(48) 
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5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 

where  P j  and  Q/  are  power  series  in  e.  It  is  seen  that  H0  is  a  function  of 
L  only  and  that  H1  is  independent  of  h.  The  notational  convention  adopted 
for  ignorable  coordinates  is  that  they  will  be  replaced  by  a  dash  in  the 
Hamiltonian's  functional  statement.  Thus  in  this  case  one  may  write 


H  =  H0  (Li  +  H1  ( L ,  G,  H,  l,  g,  -)  . 

To  eliminate  l  and  obtain  first  order  results,  one  applies  the  first  two  of 
equations  (41)  to  obtain: 


In  order  to  obtain  Sj  one  notes  that  may  be  separated  into  a  secular  part 
independent  of£  (H1  )  and  a  periodic  part  dependent  upon  l  i.e. 


f /  1  3  H* 2  \ 

L' 3  1 

-  I  -  T  +  o 

(52) 

.  • 

H  (L*.  G\  H',  g ,  -) 
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5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 


H  (L1  ,  G\  H 1 ,  l,  g, 


U4fe2 


I  6 


H 


i  2 


J_  3 _ 

2  2  G ' 2 


(53) 


_3_ 

2 


3_ 

2 


H  2 


■  I  2 


Thus  one  may  write 


and 


h;(l 

,  G',  H1, 

■»  9’  ~) 

=  «1$  (L. 

,  G  1  ,  H' ,  -,  g ,  - 

)  = 

u“fe2 

[/  1  +  3  H'2\ 

L* 3  ' 

L'6 

\  2  2  G ' 2  / 

G  13 

3S! 

-  "X 

(  i-  + 

3  H'2 

\  .  + (±__L 

H'2  \ 

3£ 

L'3 

\  2 

2  G'2 

/  6i  \  2  2 

G'2  / 

(54) 


(55) 


or 


Si  = 


P  fe. 


I  3 


/  1  3  H'2  \  f  jo  /  3  3  Hl2  \  /* 

\  J*T  T77) J  6,dt+ 


s2  de 


.(56) 


Using  equations  (47)  -  (48),  the  last  expression  becomes 

u2fe. 


-  —  +  —  — \  V  —  p  j  sin j£  +  ) 

2  2  G ' 2  ;  \  2  2  G  ’ 2  / 


2  G 


2  G' 


y-i 


aj 


sin  (2g  +  jl) 


i  =  -°° 


(57) 
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5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 

Equations  (39)  -  (40)  may  be  used  to  obtain  expressions  for  L,  G,  H,  l,  g,  h 
through  first  o^der  from  equation  (57).  The  last  of  equations  (41)  can  then 
be  applied  to  obtain  second  order  expressions  and  the  entire  process  repeated 
to  eliminate  g  and  obtain  the  final  results. 

This  is  essentially  the  procedure  followed  in  the  complete  development  of 
the  Brouwer-Lyddane  theory  3>4  except  that  there  the  significant  effects  of 
the  third,  fourth,  and  fifth  zonal  harmonic  terms  of  the  geopotential  are 
also  included.  The  results  of  the  Brouwer-Lyddane  theory  as  used  in  PULSAR 
are  summarized  below.  Define: 


a"  =  semi -major  axis  constant  (input) 
e"  =  eccentricity  constant  (input) 
l"  =  inclination  constant  (input) 
t  =  time  from  epoch 
no  -  b'a"*)h 

n  =  (1  -  i"*)h 

0  =  cosi " 

Y2  =  y  C20  a\  /a"2 

Y1  =  "4 

2  y2  n 


-C30  *3  a"-3  n-6 


u  r  _  4  _  II 

"  A  L ,  .  CL 

8  e 


"  ~^50  a"'5  n”10 

a  =  1  -  502 


> 


(58) 


3  Brouwer,  D.,  "Solution  of  the  Problem  of  Artificial  Satellite  Theory 
Without  Drag",  The  Astronomical  Journal,  Vol.  64,  No.  1274,  1959,  pp.  378-397. 

4  Lyddane,  R.  H.,  "Small  Eccentricities  or  Inclinations  in  the  Brouwer  Theory 
of  Artificial  Satellites",  The  Astronomical  Journal,  Vol.  68,  No.  8,  1963, 
pp.  555  -  558. 


ilk.  Tk.i 


i 


i 


■  -1 

.  H 

..  t 
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2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 

3=1-  lie2  -  4O04  a-1 
Y  =  1  -  302  -  804  a-1 

5  =  i  _  902  -  240 4  a-1 

X  =  1  -  502  -  166 4  a'1  . 

en  the  secular  terms  are  computed  from: 


nQ  t  | 1  +  y  y2  1  n(302  -  D  +  ~^T  2  n  |  "15  +  16n 
+  25n2  +  (30  -  96n  -  9On2)02  +  (105  +  144n  +  25n2)  64  ] 


+  ilY^ne"2  (3  -  3O02  +  3564  ]  j  +  l0"  , 


(58) 

(con't) 


n  t  i  -  —  y  '  a  +  —  y  '  2  1-35  +  24n  +  25n2  +  (90  -  192n 
V  )  2  i  32  2 


32  ’2 


-  126n2)02  +  (385  +  360n  +  45n2)04]  +  -^rY4’  [  21  -  9n2 
+  (-270  +  126n2 )62  +  (385  -  189n2)04]|  +  So"  > 


(  3  2 

h"  =  nQt  j-3y2  0+yY2' 


|(-5  +  12n  +  9n2)0  +  (-35  -  36n  -  5n2)03] 


+  —  Y  '  (5  -3n  )  6  (3  -  702)  >  +  h0" 
4  4  J 


ie  long  period  (dependent  upon  g")  terms  are  computed  from: 
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5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 


6'e  =  T6‘T7e"Vx  s1ai" K'28-^'t) 


35  Y  1 

sin2g"  — 5— e"2n2 


e"2  n2  A  siru."  sing"  +  — 


V  +  -jg- Y,.'  (4  +  3e"2)  S |sin+"  sing"  +  |3y2 


24y*  2 


3y  1  B  -  10y  "y 


V+g>  =  g"  +  i»  +  -L  ;  ~3y2  1 2  1 2  +  e"2  -  11(2  +  3e"2)02  -  40(2  +  5e"2) 

04  a-1  -  400  e"206a"2|  +  10y  1  |  2  +  e"2  -3(2+3  e"2)  02 

-  8(2  +  5e"2)  eV1  -  SOe"2  06cT2  J  +  ~~  111  e_  A  ^  -  y  J  | 


sin  2g"  + 


35  Y,  ' 


5 _  ^  3  _ji 


384  y2  ' 


n  e"  A  siru"  + 


35  Y. 
1152 


e"  302 


A  -e"  (3  +  :_"2)  siru"  +  — -J  +  2e"  3  82  siru" 


5  +  3202a-1  +  809  V2|  Los  3q"  +  I  -  +  J_  V 

1  ‘J*  (  4y2 1  siru"  64  y  1 


“  e"  (4  +  3e"2)  +  e"  siad"  (26  +  9e"2)  6  -  — 


32  Y; 


e"  92  sinT"  (4  +  3e"2  )  (3  +  1602a_1  +  400 V 2 )  +  ~  H-  siru"' 

4  Y  ’ 

2 
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,2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 


(  — — A  (3-e"2  (3  -  e" 2 ) j +  4-~T- 

\  1  +  n3 /  '  ‘  64  y. 


n  6 


e" (-32  +  81e"4) 


4  +  3e"2  +  n (4  +  9e"2) 


(63) 
(con ' t) 


siru  >  cosg 


id 


35y  5'  e"30 


h'  =  h“  + 


144  y  2 


AsinV"  +  siru"  I  5  +  320  2a  1  +  800 4oT 2 


II 


e"20 


sin2g"  cosg"  + 


12y  1  V 


1  2 


11  +  8O0V1  +  2OO04cf2  +  10y 


3  +  16G  V 1  +  400 4a” 2 |  )  sing"  cosg"  + 

|  jX  sin "U"  +  siru"  (5  +  320  V1  +  8O0V2) 


35y 


576y2 1 


5-  e" 3  g 


(64) 


e"0 


4y  'siru" 

2 


Y  *  + 

3 


5  Y  '  (4  +  3e"2)<5  +  — —  y  '  (4  +  3e"2  )(3  +  160V1  +4O04a  2) 


16 


8  s 


sinV 


+  "  J |  cosg" 

he  short  periodics  (dependent  upon  E',  f',  £")  are  computed  from: 


a  =  a"  -  a"  ~  (302-l)  + 
^,3 


a 


(1  -  e"cosE')3 


39  -  1  +  3  sinu 


2  ;n 


cos(2g"  +  2f) 


(65) 
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5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 


e  =  e"  +  5  e  + 


(  3 ':-l 


2  \ 


e"  j 3  -  e"2 (3  -  e"2)  j  +  j  3  +  e"cosf' 
1+n3  I 


(3  +  e"  cosf 1 ) >  cosf 


3(1  -  9 2 ) 


e"  +  |  3  +  e"cosf'  (3  +  e"cosf')[  cosf 


cos(2f 1  +  2g")J  -  — J—  (1-02)  3cos  ( 2g "  +  f •  )  +  COS ( 2g "  +  3f '  ) 


(66) 


e"e 


n2sirvc" 


e  +  e"  y  1  0  sinf."  sinf'  sin(2f'  +  2 g")  +  2e"Y2 ' 9sin-c" . 

(67) 


cosf1  cos(2f'  +  2 g")  +  —  Y  ‘  9  sirvt"  cos  ( 2f 1  +  2g") 


g  +  l  =  g1  +  V  + 


| -6a  (f*  -  £"  +  e"  sinf')  +  (3  -  502)  [ 3s i n ( 2f 1 


+  2g" )  +  3e"  sin  (2g"  +  f ')  +  e"  sin  (2g"  +  3f')]j 


e"  n^Y2' 


4(1  +n) 


2 ( 30 2 -  1)  (o+  1)  sinf' 


+  3(1  -  62 )  [(1  -  a)  sin  (2g"  +  f1)  +  (o  +  1/3)  sin(2g"  +  3f')]j 


(68) 


h  =  h 1  + 


2eMY  1  9  cosf'  +  —  y  '0 

2  2  2 


sin(2g"  +  2f ' )  -  e"Y  '9  sinf'  cos 


(69) 


(2f '  +  2g")  -  3y  '9  (f '  -  t"  +  e"  sinf') 
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RESULTS  OF  THE  BROUWER-  LYDDANE  THEORY  (con't) 


l6r  =  TV  sin  2g«  Sinl" 


5_  _Y 
64  y. 


5  3 


sin\" (4  +  9e2")  6  /  •  cos  g"  + 


35  Ye 


384  Y 


—  n  e 


3e"2X 


siru"  cos3g" - -y  'n3  <2  (302  -1)  (a  +  1)  sinf1  +  3(1  -  02) 

4  2  I 


1  -  a)  sin(2g"  +  f')  +  (0+  1/3)  sin(2g"  +  3f ' ) 


i  = 


n 


1  -  e"cosE' 


1  -  e"cosE' 


:centric  anomaly  E1  is  obtained  from  a  Newton-Raphson  iteration  upon 
3pler  equation  (refer  to  SECTION  9) 


(70) 


(71) 


>  J 


E'  -  e"  si nE '  =  i" 


(72) 


ie  true  anomaly  f'  is  found  from 


sinf' 


n  s i nE ' 


1  -  e"cosE' 


:osf 1 


cosE'  -  e" 

1  -  e"cosE 1 


(73) 


inal  osculating  values  for  a,  i,  and  h  are  computed  from  equations  (65), 
snd  (69),  respectively.  Equations  (66),  (68),  (70)  and  (72)  are  used 
Iculate  final  osculating  values  for  l,  g.  and  e  from  the  following 
ions: 
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A  = 

e  cos£ 

"  -  e6£  sin£" 

(74) 

B  = 

e  sin£ 

"  +  e<5£  cos£"  , 

(75) 

l  = 

tan"1 

(B/A) 

(76) 

9  = 

U  +  9 

)  -  £ 

(77) 

and 

e  =  (A2  +  B2  \2  .  (78) 


The  inertial  position  and  velocity  are  obtainable  from  the  osculating 
Brouwer  elements  through  the  following  relationships: 


where 


cos h  cosg  -  sinfi  cos<  sing 


cos h  sing  -  sinh  cos-c  cosg 


0 


cosg  sinfr  +  cosh  cost  sing 


cos h  cos l  cosg  -  sinfr  sing 


(81) 
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5.2.2  RESULTS  OF  THE  BROUWER-LYDDANE  THEORY  (con't) 


It  should  be  mentioned  in  closing  that  the  input  Brouwer  orbital  elements 
are  obtained  from  NAVSPASUR  and  that  the  semi-major  axis  thus  provided  is 
the  Kaula  semi-major  axis  expressed  in  earth  radii.  The  Kaula  semi-major 
axis  must  be  converted  from  units  of  earth  radii  to  linear  units  using 
the  following  relation: 


a"  aK  ae 


1  + 


2X  \  K 


1  -X 


where  &K  is  the  Kaula  semi -major  axis  obtained  from  NAVSPASUR  and 


SECTION  6 


VARIATIONAL  EQUATIONS 

If  an  approximate  or  initial  orbit  is  available  for  a  satellite,  further 
orbit  improvement  involves  the  application  of  differential  correction  methods. 
Such  methods  are  used  to  improve  orbital  parameters  in  the  point  editing  and 
cross-pass  editing  processes  and  will  be  discussed  in  detail  in  SECTION  7. 
Fundamental  to  the  parameter  improvement  process  is  the  evaluation  of  partial 
derivatives  of  r  and  r  with  respect  to  those  dynamic  parameters  being  im¬ 
proved,  e.g.  orbital  elements,  drag  coefficient,  etc.  In  PULSAR  values  for 
these  partial  derivatives  are  obtained  from  the  numerical  integration  of  the 
variational  equations. 

In  the  following  subsections  the  form  of  the  variational  equations  used  by 
PULSAR  is  developed.  Also  included  in  the  last  subsection  is  a  discussion 
of  the  state  transition  matrix  which  has  important  applications  in  the  point 
editing  process  and  in  the  propagation  of  the  state  improvements  (SECTION  7). 

6.1  DEVELOPMENT  OF  THE  VARIATIONAL  EQUATIONS  AND  ASSOCIATED  ANALYTICS 

A  mathematically  rigorous  method  of  obtaining  the  partial  derivatives  of  r 
and  r  with  respect  to  those  dynamic  parameters  which  are  being  improved  in 
the  differential  correction  process  is  through  the  use  of  the  variational 
equations.  In  PULSAR  these  equations  are  integrated  numerically  along  with 
the  equations  of  motion  and  are  stored  along  with  the  ephemeris  on  the  trajec¬ 
tory  file.  The  partial  derivatives  obtained  in  this  way  include  all  of  the 
secular  and  periodic  variations  introduced  by  the  satellite's  force  environment. 
The  variational  equations  and  closed  form  analytic  expressions  for  coefficients 
appearing  in  the  equations  are  developed  below. 

6.1.1  THE  VARIATIONAL  EQUATIONS 

The  variational  equations  are  a  system  of  differential  equations  which  can  in 
general  be  derived  from  the  acceleration  function  for  the  satellite.  If  p 
is  an  orbital  or  dynamic  model  parameter  which  is  to  be  improved  via  the 
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6.1.1  THE  VARIATIONAL  EQUATIONS  (con't) 

differential  correction  process  and  is  the  parameter  at  epoch  tfl,  then 
the  acceleration  function  can  be  differentiated  with  respect  to  to  obtain 


3r 

|  P4- 

3  r 

da  dr 

-f 

3  a 

+ 

d% 

9r 

3r  3 Pk0 

» 

where 


and 
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6.1.1  THE  VARIATIONAL  EQUATIONS  (con't) 


In  the  above  expressions  a  represents  the  sum  of  accelerations  acting  upon 
the  satellite.  For  the  PULSAR  system  a  is  given  by  the  right  hand  side  of 
equation  (1)  of  SECTION  5. 


If  one  lets 


3r 

*Pk0 


(6) 


represent  the  solution  to  equation  (1),  then  it  may  be  recast  into  the  form 
of  a  system  of  inhomogeneous  second  order  1  inear  differential  equations: 


where 


9a 

and  G,  =  - 

dPk0 


(7) 


(8) 


Note  that  if  the  acceleration  function  does  not  depend  explicitly  upon  the 

parameter  then  equation  (7)  reduces  to  a  system  of  linear  homogeneous 

second  order  equations  since  G  =  0. 

1  ~ 
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6.1.1  THE  VARIATIONAL  EQUATIONS  (con't) 

The  only  dynamic  parameters  improved  by  the  differential  correction  process 
in  PULSAR  are  the  set  of  osculating  orbital  elements  (or  the  associated 
inertial  position  and  velocity  components)  and  the  drag  factor  D  =  CDA/W  , 
i  .e. 


where  a,  e,  i,  ui,  M,  and  ft  are  the  Keplerian  orbital  elements.  Consideration 
of  the  force  model  used  by  PULSAR  (SECTION  5)  provides  the  following  general 
expressions  for  the  coefficient  matrices  Gj  and  G, ,  as  well  as  the  forcing 
term  G  : 

^  3a E  9(1  D,  9<z  moon  9a  sun  ,  (10) 

=  ^  ?  If 


G 


and 


9 


_9f_0 

ao 

for 

Pfe  = 

D  =  CdA/W 

0 

for 

Pfe  = 

element  set 

(ID 


(12) 
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6.1.2  CLOSED  FORM  EXPRESSIONS  FOR  G, 


G2 ,  and  G 3 


The  coefficient  matrices  G2  and  G2,  as  well  as  the  forcing  term  G3,  can 
be  expressed  in  closed  form  from  consideration  of  the  explicit  representa¬ 
tions  of  the  accelerations  contributing  to  the  force  model  (SECTION  5). 

In  the  subsections  below,  closed  form  expressions  are  obtained  for  each  of 
the  terms  on  the  right  hand  sides  of  equations  (10)  -  (12)  above.  These  re¬ 
sults  can  be  summed  to  give  closed  form  expressions  for  the  G  matrices. 


da  F 


6. 1.2.1  Position  Derivatives  of  the  Geopotential  Acceleration: 


9r 


Let  Vj  and  Vebe  gradient  operators  in  the  inertial  and  earth-fixed  reference 
frames,  respectively.  These  operators  are  related  by  the  following  trans¬ 
formation: 


L-l 


VT  =  (B  C  D)'  -  V 
1  ^  ~  ^  “ 


where  BCD  are  the  transformations  described  in  SECTION  2  and  the  gradient 
operators  expressed  in  matrix  notation  have  the  form 


/£\ 


_3 

9y 


i 


\ 


3zy/ 


,  iy 


I,  e), 


From  equation  (2)  one  sees  that 
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6. 1.2.1  Position  Derivatives  of  the  Geopotential  Acceleration-  3r 


where  the  superscript  T  means  matrix  transposition.  Applying  equation  (8) 
of  SECTION  5  to  the  last  expression  gives 


oo  n 


■  <“££)“  Z  E  K 


n=0  m=0 

n^l 


v T  ( v  y  m ) 

e  v  e  un  'x 


v  T  (v  u") 
e  v  e  un  'y 


V  (^e  Um)z 

e  e  n  z 


V  T  (v  V  m  ) 
e  v  e  v  n  'x 


v  T  (v  V  n  ) 

e  '  e  v  n  'y 


V  <?eVnra>z 


(B  C  D)‘ 


Using  equations  (13)  -  (14)  from  SECTION  5,  it  is  seen  that 


V  1  (V  U  ) 
e  v  e  un  'x 


V  T  (V  U™  ) 
e  v  e  un  'y 


V„T  (v„  II m). 


m  +  T  m-1  ->  t  m+1 

A  V  U  -  V  u 
Hn  e  Vl  e  un+1 


m  ->T  m-1  +  t  m+1 

-A_  V  V  -  V  V  n 

n  e  n+1  e  n+1 
-2(n-m+l)  vJ  U™. 
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6. 1.2.1  Position  Derivatives  of  the  Geopotential  Acceleration:  r  (con't) 


V  T  ( v  V  m  ) 
e  '  e  vn  ;x 


v  T  (v  V  n  ) 

e  v  e  ¥n  'y 


V  (vevnm)z 


m  ->  t  m-1  y  m+1 

A  V  V  -  V  y 
Mn  e  v  n+!  e  vn+i 


m  -+T  m-1 


An  VeUn+1  +  Vll 

\  -2(n-m+l)  V  '  V 
\  e  n+1 


where 


l\n  =  (n-m+l)(n-m  +  2)  . 

The  gradients  of  the  and  terms  in  the  right  hand  sides  of  equatii 

■i  -L 

(17)  -  (18)  above  are  themselves  computed  from  equations  (13)  -  (14)  of 
SECTION  5. 

9ct  | 

6. 1.2. 2  Position  Derivatives  of  the  Atmospheric  Drag  Acceleration:  9r 
From  equation  (16)  SECTION  5,  it  is  seen  that 


-TD9 


9P  -»■  9v  ->  9v 

— — -  vv  +  p  — rr-  V  +  pu  - 

9r  3r  9r 


Consider  the  first  partial  derivative  on  the  right  hand  side  of  the  last 
equation: 

9p  3p  9h 

9r  9h  9r 


From  equation  (19)  of  SECTION  5,  one  finds  that 
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6. 1.2. 2  Position  Derivatives  of  the  Atmospheric  Drag  Acceleration:  —+r  (con't) 

3r 


3p  1  2C3^  +  C4 

3h  1  2  (C-fi2  +  zh  -  r)** 


Equation  (20)  of  SECTION  5  may  be  rewritten  as 


r  1  - 


2  (  £2  \  2  ** 

r  +  (  1  -  e2  )  X  3 


where  is  the  semi  major  axis  of  the  reference  ellipsoid,  e  is  the  earth's 
eccentricity  and  r  is  the  magnitude  of  the  satellite  inertial  position  vector. 
Note  that  the  notation  (xi,  x2,  x3)  =  (x,  y,  z)  has  been  adopted. 

i.  u 

Taking  the  partial  derivative  of  h  with  respect  to  the  l  position  component 
and  rearranging  results  gives 


(r  -  h)  3  (r  -  h)3  e2 

- —  +  — - T  6 


(1  -  e2) 


where  6  is  the  Kronecker  delta.  This  result  may  be  expressed  more  generally 


*11 


(r  -  h)3  vr  ,  (r  -  h)3  e2  ^ 

h  + - r  +  -  r  e  I  e 

a  e  a  U  -  e  ) 

C  J  e 


-96- 


TR  82-391 


da  q 


6. 1.2. 2  Position  Derivatives  of  the  Atmospheric  Drag  Acceleration:  3r  (con‘t) 


where  e^  is  the  unit  vector  along  the  inertial  x3  axis. 

Consider  now  the  second  partial  derivative  on  the  right  hand  side  of  equation 
(20).  Since 


/-*■  T  -+\k 

-  (V  V) 


then  by  taking  the  partial  derivative  of  equation  (26)  with  respect  to  the 
^th  position  component,  one  obtains 


(26) 


9v 

3X; 


V 

V 


3x, 


(27) 


or  more  generally 


3  v 

9  r 


3  r 


(28) 


3v 


An  expression  for  3r  is  now  needed.  Once  this  is  obtained,  both  equation 
(28)  and  the  third  partial  derivative  on  the  right  hand  side  of  equation  (20) 
will  be  known  in  closed  form. 


The  velocity  of  the  satellite  relative  to  the  rotating  atmosphere  is  given 

t  h 

in  Section  5  by  equations  (17)  -  (18).  Thus  the  relative  velocity's  i 
component  is  given  by 
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3a  D 


.2.2  Position  Derivatives  of  the  Atmospheric  Drag  Acceleration: 3r  (con't) 


v  . 
A. 


=  -  (w  X  r)^  =  L  -  e ;;t.  CO;  xt,  U  =  1,2,3)  , 


-ijk  Afe 


y-i 

fe=i 


re  e..,  is  the  Levi-Civita  density  defined  by 


(29) 


<!*■ 


0  if  any  index  is  equal  to  any  other  index 
+1  if  x. ,  j,  k  form  an  even  permutation  of  1,  2,  3 
-1  if  -c,  j,  fc  form  an  odd  permutation  of  1,  2,  3  ■ 


(30) 


ing  the  partial  derivative  of  equation  (29)  with  respect  to  the  l 
ition  component  gives 


th 


9v . 


)X£ 


3X£ 


E 

y- 1 

fe-i 


£  .  0)  . 
J 


3xfe 


3 

E 

y=l 

fe-1 


%-k  w j  6kZ 


(31) 


9  x 


£ 


E 

y=i 


=  -  >  e  -  n  to  - 

j 


(32) 


?se  results  may  be  collected  to  give  the  desired  closed  form  result: 
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6. 1.2.2  Position  Derivatives  of  the  Atmospheric  Drag  Acceleration:  9r  (con't) 


4°)  =  It  .  >  dJ/HN  J_  „.x.  k  *  ^ 

3r  yL  3Xy  2  (  \3A  /  r2  <•  J  ag2 


(r-k)  e2  p  3v  h 

+  — ; -  6-  +  —  V.  )  v ,  11  +  pv 

a  2(l-e2)  J 3  v  *  L-,  k  9X; 


ae2(l-e2)  J3 


v  ^  ,  vfe  3X 


where  the  position  derivatives  of  velocity  are  given  by  equation  (32) 


6. 1.2. 3  Position  Derivatives  of  the  Luni-Solar  Gravitational  Accelerations 


ar  ,  \ z  =  sun,  moonj 


The  luni-solar  gravitational  acceleration  is  given  by  equation  (26)  in 

J.  L_ 

SECTION  5.  Taking  the  partial  derivative  if  the  acceleration  component 
with  respect  to  the  position  component  gives 


d(«ih 


-  \.)  *t.  ) 

- - - — T rr  +  r - T-T7T  /  >  U  =  sun,  moon).  (34) 
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7.1.1  APPLICATION  TO  THE  POINT  EDITING  PROCESS  (con’t) 

pass  normal  equations  for  the  parameter  corrections  are  outined  in  SECTION  9. 

7 . 1 . 1 . 1  Computed  Observations  and  Observation  Partial  Derivatives 

It  can  be  seen  from  equations  (14)  -  (15)  that  in  order  to  form  the  8  and 
E  matrices,  one  must  not  only  be  able  to  provide  computed  values  for  the  ob¬ 
servations  to  be  processed,  but  one  must  also  be  able  to  form  the  matrix  of 
data  partial  derivatives  with  respect  to  the  parameters  to  be  improved.  Fur¬ 
thermore,  it  nay  be  necessary  to  correct  the  computed  observations  for  effects 
introduced  by  variable  attitude  and  displacement  of  the  transmitting  antennae  from 
the  center  of  mass.  The  computed  observations  must  also  be  corrected  for  tropospher¬ 
ic  refraction ,  antenna  rotation  and  general  relativistic  effects.  Details  concerning 
the  observation  corrections  can  be  found  in  SECTION  4.  The  time  associated 
with  each  observation  has  previously  been  corrected  for  transmission  delay  and 
clock  errors  before  the  observations  are  edited  by  the  point  editor. 

The  data  processed  by  the  point  editor  are  range  difference  observations.  A 

t  h 

computed  value  for  the  -t  n  range  difference  point  is  obtained  from 

C.  =  V(Pj,  t.)  -  V(p.,  t.  .)  +  6AR-,  (16) 

-v  J  A,  j  I  A, 


where 


V{pj,  t) 


fb  -  V 


2f  *b  ^  "  M2  +  ^  +  CR)AfR 


(17) 


and  o,‘  R .  is  the  correction  for  antenna  rotation  which  is  discussed  in  SECTION 

L 

4.  It  is  considered  here  as  a  correction  to  the  observed  range  difference,  but 
has  been  added  to  the  computed  range  difference  for  convenience.  Consequentl 
its  partial  derivatives  will  vanish  when  forming  the  a  matrix  | equation  (19) 


Since  ] ,'Tr^ j  is  small,  one  may  Taylor  expand  the  exact  expression  for  |p 
given  by 
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(14) 

(15) 

.1  APPLICATION  TO  THE  POINT  EDITING  PROCESS 

■  PULSAR  applications  a  pass  of  tracking  data  for  a  given  satellite  is 
ined  to  be  the  set  of  observations  made  by  a  s;ng!e  tracking  station  during 
;  of  its  satellite-station  inview  periods  for  the  satellite.  Since  all  ob- 
•vations  in  a  given  pass  are  not  necessarily  of  equivalent  accuracy,  it  is 
:essary  to  identify  observations  which  are  inconsistent  with  other  observa- 
>ns  in  the  same  pass  so  they  will  not  corrupt  results  obtained  by  any  post- 
_SAR  processing  of  the  observation  data.  This  identification  procedure  is 
'formed  by  PULSAR'S  point  editing  process. 

2  point  editor  processes  observations  on  a  pass-by-pass  basis  by  fitting 
2  satellite's  orbit  to  a  pass  of  data  using  the  method  of  weighted  least 
jares.  The  resulting  fitting  parameter  improvements  Ap  are  used  to  adjust 
2  observation  residuals  0_-_C .  The  quality  of  the  adjusted  residuals  is 
sessed  and  if  the  quality  of  certain  of  them  is  not  acceptable,  those 
nervations  are  tagged  as  outliersand  they  are  not  used  in  the  formation 
the  normal  equations  for  the  next  iteration  of  the  associated  differential 
r’rection  process.  However,  the  quality  of  the  residuals  associated  with 
2  tagged  observations  are  assessed  during  each  iteration  of  the  correction 
acess  and  if  their  quality  is  found  to  be  acceptable  during  these  iterations 
2  data  points  are  untagged  and  used  in  the  formation  of  the  normal  equations 
subsequent  iterations  of  the  differential  correction  process. 

this  subsection  only  the  underlying  mathematics  used  in  forming  the  pass 
rmal  equations  and  the  associated  differential  correction  process  are  dis¬ 
used.  The  statistical  methods  used  to  assess  the  quality  of  the  observations 
d  the  criteria  used  to  define  the  iterative  convergence  conditions  are 
dressed  in  SECTION  8.  The  numerical  techniques  used  in  the  solution  of  the 
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7.1  THE  METHOD  OF  WEIGHTED  LEAST  SQUARES  (con't) 


confidence  levels  between  the  various  observations,  as  well  as  possible 
inter-relationships  between  the  observation  errors).  This  statement  is 
equivalent  to  a  weighted  least  squares  minimization  procedure,  i.e.  find 
the  value  of  <A p>  which  will  minimize  the  functional  J  given  by 

J(<Ap>)  =  |(0  -  C)  -  a<Ap> J  ^  W  j  (0  -  C)  -  a<Ap>j ,  (? 

where  W  is  an  M  X  M  symmetric,  positive  definite  weighting  matrix  which 
is  assumed  to  be  known.  Taking  the  variation  of  J  with  respect  to  <Ap>  gives 

SJ  «  -  6<Ap>T  aT  W  [(0  -  C)  -a<Ap>|  -  £(0  -  C)  -a<Ap>J  T  W  a  6  <Ap> 

=  26<Ap>T  aT  W  [(0  -  C)  -  a<Ap>  J  ,  (] 

provided  that  W  =  WT.  If  the  above  equation  is  to  be  an  extremum,  it  must 
be  zero.  Since  the  variation  6  <Ap>  is  completely  arbitrary,  one  has 

J  U  [  (0  -  C)  -  a<Ap>  =  0  ,  (] 


(a  W  a)  <Ap> 


a  W  (0  -  C). 


This  equation  is  referred  to  as  the  normal  equation  for  the  weighted  least 
squares  problem.  In  subsequent  sections  the  best  estimate  notation  will  be 
suppressed  and  the  normal  equations  will  be  written  as 


where  it  is  understood  that  Ap  is  the  "best  estimate"  in  the  weighted  least 
squares  sense  and  8  and  E  are  given  by 
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7.1  THE  METHOD  OF  WEIGHTED  LEAST  SQUARES  (con't) 

1L 

Retaining  only  the  linear  terms,  one  then  has  for  the  observation 


(0  -  C) • 


N 

£ 

r  i 


9 V(Pj,  tj) 
9  Pi 


Pj  =  P°i 


For  all  M  observations  the  following  matrix  equation  may  be  written: 
(0  -  C)  =  a  Ap  +  e  , 


(6) 


(7) 


where  (0_^C)  is  the  M  X  1  residual  matrix,  a  is  the  M  X  N  matrix  of  partial 
derivatives  of  the  observations  with  respect  to  the  unknown  parameters  py, 

Ap  is  the  N  X  1  matrix  of  parameter  corrections,  and  e  is  the  M  X  1  matrix 
of  random  errors. 

As  a  result  of  the  last  expression,  the  orbit  estimation  problem  may  be 
restated  as:  given  the  matrix  of  observation  residuals  (O^^C)  and  the 
matrix  a,  find  the  "best  estimate"  of  Ap,  i.e.  <Ap>.  Once  <Ap>  is  known, 

«v  /V 

the  "best  estimate"  of  the  unknown  parameters  py,  i.e.  <py>  is  obtained  from 

<Pj>  =  P°j  +  <Apy>  ,  j  =  1,  2,  ...,  N.  (8) 

It  should  be  noted  that  since  the  estimation  problem  has  been  linearized 
(i.e.  only  the  n  =  1  term  retained  in  equation  (5)), equation  (7)  may  have 
to  be  solved  many  times  (in  the  iterative  sense)  to  obtain  a  solution  to  the 
nonlinear  problem.  The  method  of  nonlinear  least  square  curve  fitting  is 
called  differential  correction. 


The  "best  estimate"  in  the  least  squares  sense  will  be  that  particular  set 

of  values  <Ap>  which  minimizes  the  weighted  sum  of  squares  of  the  elements 

of  the  residual  matrix  (weighting  is  used  to  account  for  the  difference  in 

» 
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7.1  THE  METHOD  OF  WEIGHTED  LEAST  SQUARES  (con’t) 

°i  -  ^  •  ‘2> 

J.  L. 

where  the  quantity  represents  the  error  associated  with  the  -t  observa¬ 
tion.  It  should  be  noted  that  it  is  assumed  that  no  errors  exist  in  the 
time  at  which  a  particular  observation  is  being  made. 

In  general  V(pj,  t)  is  nonlinear  in  the  parameters  pj .  Even  though  partic¬ 
ular  nonlinear  functional  forms  could  be  applied  to  specific  cases,  it  is 
desirable  to  obtain  a  method  which  is  applicable  in  the  general  case.  To 
facilitate  this  suppose  it  is  known  that  the  true  values  of  the  pj  lie  close 
to  some  prescribed  set  p°j .  Then  the  values 

C.  =  V(p°j,  t^)  (3) 

can  be  computed.  The  observation  residual  for  the  observation  can  be 
formed  by  differencing  equations  (2)  -  (3)  to  give  the  observation  residual 

(0  -  C)  .  =  0  .  -  C  .  =  V(pj,  t  )  -  V(p°j  ,  t  )  +  e  .  .  (4) 

4,  4.  -v  A,  A, 


It  should  be  noted  that  the  observation  residual  is  still  corrupted  by  the 
errors  arising  from  the  observation  technique. 


Expanding  the  right  hand  side  of  the  last  equation  in  a  Taylor  series  about 
Pj  gives 

CO  N  -,i 

(O-C),  .  ptf.  V  ♦£(±,2 

1=1 


n=l  n!  ”  9Pi 


PJ  *  PJ 


(t-y  -  p )■)"  +  e€-  -  viPj,  t.j. 
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SECTION  7 

STATE  ESTIMATION  AND  PARAMETER  IMPROVEMENT 

Since  a  very  large  number  of  observational  measurements  of  an  earth  satel¬ 
lite's  position  and  velocity  are  usually  available,  the  satellite  orbital 
state  estimation  problem  is  generally  highly  over-determined.  This  remains 
true  even  when  various  modeling  parameters  (e.g.  drag  coefficient,  station 
parameters,  etc.)  are  introduced  into  the  problem.  As  a  result,  the  state 
estimation  process  is  very  effective,  since  it  can  be  supported  by  least 
squares  estimation  methods  which  also  permit  a  statistical  interpretation 
of  the  orbit. 

In  the  following  subsections  the  method  of  weighted  least  squares  is  discussed 
and  the  associated  normal  equations  are  developed.  These  results  are  then 
applied  in  a  description  of  the  processing  which  is  utilized  by  both  the  point 
editor  and  cross  pass  editor  of  PULSAR.  Also  included  is  a  brief  discussion 
concerning  the  state  improvement  propagation  technique  used  in  PULSAR. 

7.1  THE  METHOD  OF  WEIGHTED  LEAST  SQUARES 

The  fundamental  orbit  estimation  process  is  that  of  solving  for  a  set  of  obser¬ 
vation  model  parameters  X/  in  such  a  manner  that  the  constrained,  weighted  sum 
of  the  squares  of  the  differences  between  the  observation  model  and  the  observed 
trajectory  is  minimized.  To  initiate  the  development  of  the  mathematics  related 
to  the  orbit  estimation  problem  consider  an  arbitrary  curve  in  space  which  is 
described  by 

ij(t)  =  V(p p2»  ...,  pN,  t)  =  P(pj,  t)  (1) 

where  the  pj  (j  =  1 , 2, . . . ,N)are  as  yet  unspecified  parameters.  Assume  at  the 
times  ti  {l  =  1 ,  2,  . . . ,  M;  M  >  N)  observations  are  made  of  y{t-).  In  general 
the  true  value  of  yit^)  will  be  corrupted  by  errors  arising  from  the  obser¬ 
vation  techniques,  so  that  the  observation  at  time  t.  and  denoted  by  0  .  is 

■4-  4. 
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6.2  THE  STATE  TRANSITION  MATRIX  (con't) 

The  matrix  <T‘(t  )  $  (t  )  is  referred  to  as  the  state  transition  matrix, 

since  it  enables  the  variational  equation  solutions  referenced  to  one  epoch, 
say  t0,  to  be  transformed  to  another  reference  epoch  tx. 


-107- 


TR  82-391 


6.2  THE  STATE  TRANSITION  MATRIX  (con't) 

Pfe( t o )  re^ers  to  the  orbit  elements  at  epoch  t0.  Thus  one  may  write  for  the 
general  solution 


¥  (t)  =  $  (t)  K 

-PfcftJ  ~  ~ 


(61) 


where 


¥  (t) 

~t 


0 


(l  <t> 

\  p,(t.) 


¥  (t)  ¥  (t)  ¥  (t)  ¥  (t)  ¥  (t)  \ 

'P2(t0)  'P3(t0>  \(t.)  'Ps<t0) 


Solving  for  K  at  t  =  t}  gives 


K  »  4  (t  )  ¥  (tj 

“t0  "PfeltJ 


(63) 


so  that  equation  (61)  becomes 


¥  (t)  =  $  (t)  Mt  )  ¥  (t,) 

“Pfe^J  ~t0  ~Pfe(t1) 


(64) 


This  expression  may  be  generalized  so  that  the  solutions  for  all  orbit 
parameters  p^  can  be  transformed  simultaneously  by  using  the  form  of  equa¬ 
tion  (62)  in  the  last  expression  to  obtain 


*  (t) 


$  (t) 

A. 


(t  ) 


! 


Equations  (43) 


(54)  are  used  to  evaluate  ¥  (tt)  and  $  (t  ). 

^felti)  \ 


(65) 
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$  D 

6.1.3  INITIAL  CONDITION  PARTIAL  DERIVATIVES:  pk  (con't) 


6.2  THE  STATE  TRANSITION  MATRIX 

It  is  necessary  to  develop  a  set  of  orbit  parameter  transformations  which 
can  be  used  to  change  the  epoch  of  the  Pfe^  of  the  variational  equation 
solutions  from  the  trajectory  epoch  t0  to  the  new  epoch  tj.  To  facilitate 
this  one  can  write  equation  (7)  in  a  more  convenient  form.  If  one  lets 


Y  (t) 

~Pfc 


\ 

then  equation  (7)  becomes 


$pfe(t) 


V‘> 


r(t)  y  (t)  +  A ( t ) 
-Pfe 


where 


0 

I  \ 

and 

A(t) 

(  0  ' 

Si 

G’  1 

G  3 

1 

When  Pfe  is  an  orbit  parameter,  the  variational  equation  is  homogeneous  in 
form  and  its  general  solution  in  terms  of  orbit  elements  at  some  epoch  tx,  d 
noted  by  Y  (t)  ,  is  a  linear  combination  of  the  Y  (t)  ,  where  the  subscript 
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1 


I 


*1 


6.1.3  INITIAL  CONDITION  PARTIAL  DERIVATIVES: 


(con1 t) 


^M+cj 


and 


t 


where  E  is  the  eccentric  anomaly  and 


Y  =  (1  -  e2)h 

P  =  a( 1  -  e2) 


(53) 


(54) 


(55) 


It  should  be  noted  that  if  position  and  velocity  are  used  instead  of  orbital 
elements,  then  the  initial  conditions  are  given  by 


where  6^a  is  the  Kronecker  delta,  i  =  x,  y,  z  and  j  =  x,  y,  z. 


(56) 


To  facilitate  integration  of  the  inhomogeneous  form  of  the  variational  equa¬ 
tions,  the  following  conditions  are  imposed  at  the  integration  epoch  for  p^  =  D: 
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$  $ 

6.1.3  INITIAL  CONDITION  PARTIAL  DERIVATIVES:  Pfc,  Pfc  (con't) 


+  w 


resinu) 


Y2"  ( — )  {~r  +  e  cos  E)  cos  +  E)  ’  e  cos 


e"  cos 


os2  E  +  y\ 
1  +  Y  /. 


-►  1  .  ,  ^  /esinE 

r  +  — —  sin(w+E)  -  e  cosu)( - — 

y2  \i +  y 


;  sin  e\|  4 
1  +  Y  /  r  ’ 


'ecosw 


n  /  a 
"y5"  \  r 


+  e  cos  E 


\  /e2  cos2  E  +y\  1 

J  sin  (w+E)  +  e  si  noil — - - - - )  I  . 


r  +  — y  cos  (oj  +  E)  +  e  sinai 


/e  sin  e\  i 

i  nai  ( - )  r  , 

\  1  +  Y  /. 


z  sin  fi 


-z  cos  n 


y  cos  Q  -  x  sin  Q 
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6.1.3  INITIAL  CONDITION  PARTIAL  DERIVATIVES:  Pfe  ,  Pfe  . 

In  order  to  integrate  the  variational  equations,  one  must  have  initial 
conditions  for  the and  vectors.  These  are  obtained  for  the  orbital 
elements  by  evaluating  the  following  closed  form  expressions  at  the  inte¬ 
gration  epoch.  Since  the  derivatives  of  these  relations  have  been  developed 
elsewhere1,  only  the  results  will  be  given  here: 


a 


r 

~a 


(43) 


resinto 


e  sin  E 

sin(oo+  E)  +  esinw  -  ecosto  — - 

1  +  Y 


r  + 


1 

n 


r  r  (e  cosE  +  Y) 

-  (1  +  —  )  cos  (o>  +  E)  +  —  e  costo - - - 

p  P  1  +  Y 


r,  (44) 


1 

'  ecosto  Y2 


cos(u)  +  e)  +  e  cos  to  +  e  sin  to 


e  sin  E 
1  +  Y 


r  + 


/  r  \  ,  .  r  (e  cos 

(1  +  — Ism  (to  +  E) - e  smto - 

V  p  /  p  1  + 


(e  cos  E  +  y) 


Y 


r  ,  (45) 


j  l  sin  ft 
-I  cos  ft 
y  y  cos  ft  -  x  sin  ft 


(46) 


1  Hubbard,  E.  C.,  Orbit  Elements  Determinate  At  Zero  Eccentricity,  NWL  Tech¬ 
nical  Memorandum  No.  K  -  26/63,  Dahlgren,  Virginia,  19^3. 
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6. 1.2.3  Velocity  Derivatives  of  the  Atmospheric  Drag  Acceleration:  (con't) 


Using  equation  (26),  one  sees  that 


I  r\  ' 

V  OV 

v  a? 


where  v  is  given  by  equation  (29).  Taking  the  partial  derivative  of  equation 
(29)  with  respect  to  the  c  velocity  component  gives 


_ i-  _  v 

*1  7T 


h,k  “y 


Collecting  these  results  gives  the  desired  closed  form  result: 


3£D\  _ 


1  /  P  \  .  2  .  1 

—  Dq  I  —  )  v .  v  .  +  v  6-.  I 

2  9  \  u  )  <  i  V  j 


6. 1.2. 5  Partial  Derivative  of  Atmospheric  Drag  Acceleration  with  Respect  to 

Bln 

the  Drag  Coefficient:  9D 

To  obtain  an  expression  for  this  forcing  term,  one  merely  differentiates 
equation  (16)  of  SECTION  5  with  respect  to  the  drag  coefficient  D  to  obtain 


-  —  g  p  v  v 

2 


where  u  is  given  in  SECTION  5  by  equations  (17)  -  (18) 
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6. 1.2.3  Position  Derivatives  of  the  Luni~Solar  Gravitational  Accelerations 


9r  ,  (£  =  sun,  moon)  (con't) 


one  finds  that  equation  (34)  becomes 


9  -up  (  ) 

-  =  K  -  3  ( X .  -  X.  )  (X.  -  X»  )  ,  (£=sun,  moon)  (36) 

9  Xj  n5  \  1  ^  S 


where  is  the  Kronecker  delta  and  l,  j  =  1,  2,  3.  Putting  this  result  in 
the  form  of  equation  (2)  gives: 


-  3<x,  -  x<t  -3(x; hx^x^  ) 


£  -VI 


-*\-h  n\-\  )  \  >2  -3(VX*  )(x,-x*  > 


-3(X  -  Xt  )(X  -  Xt  )  -3 ( X  -  xe  )(X  -  X(  )  n!  — 3 ( X  -  X;  ) 

3  1  3  2  &  3 


where  J,  =  sun,  moon. 


6. 1.2. 4  Velocity  Derivatives  of  the  Atmospheric  Drag  Acceleration:  9r 


Since  the  atmospheric  density  is  independent  of  velocity,  the  velocity  deriva 
tives  of  equation  (16)  SECTION  5  are 
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7. 1.1.1  Computed  Observations  and  Observation  Partial  Derivatives  (con't) 
to  obtain  the  nearly  exact  result  used  to  form  the  computed  observations: 


IpI  “  l^(te)  ~rs(tR)  I  +  ArA(te)  • 


r(te)  -  rs(tR) 
|r(tg)  -  r$(tR)| 


Ap 

GR  • 


(18) 


In  equation  (17), 


CR 

AfR 


and 


106  , 

frequency  bias, 
frequency  drift  bias, 
refraction  bias, 

tropospheric  refraction  correction  (see  SECTION  4), 
time  of  first  observation  of  pass. 


C  =  speed  of  light  in  vacuo, 

where 

t  =  time  of  signal  emmission  from  satellite, 

tR  =  time  of  signal  reception  at  observing  station, 

r$(t)  =  inertial  station  position  at  time  t, 

r(t)  =  inertial  position  of  satellite  center  of  mass  at  time  t, 
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inertial  correction  vector  for  antenna-  offset  and  satellite 
attitude  at  time  t  (see  SECTION  4), 


7. 1.1.1  Compute 
ArA(t)  = 
and 

ApGR  =  general  relativistic  correction  to  range  (see  SECTION  4). 

For  range  difference  observations  the  elements  of  the  a  matrix  are  obtained  by 
evaluating  the  expressions  for  the  partial  derivatives  of  the  computed  obser¬ 
vations  with  respect  to  the  parameters  to  be  improved.  Taking  the  partial 
derivative  of  equation  (16)  with  respect  to  the  parameter  Pfe  gives 


dCl  dViPj>  xi)  (P/>  Xi- 1^  3  ( 6AR^_) 

9Pfe  9Pfe  9  Pk  9  pk 


(19) 


where  V(pj,  t)  is  given  by  equation  (17).  If  Pfe  is  a  dynamic  parameter  (i.e. 
an  orbital  element  or  drag  parameter),  one  may  write  for  the  observation  partial 
derivative  at  time  t 


dD(pj,  t)  9 V[p-,  t)  3r 

9  Pfe  9r  9 pk 

or 

9 Vipj,  t)  _ 

- — -  =  VV(Pj,  t)  •  Vpfe(t)  +  VV(p-,  t)  .  Ypfe(t)  ,  (21) 

where  'Fpfe(t)  and  H'pfe(t)  are  the  associated  solutions  to  the  variational 
equations  (see  SECTION  6).  By  inspection  of  equation  (17)  it  is  seen  that 
V(pj,  t)  is  independent  of  velocity,  so  that  one  may  write 


9 V(pj,  t)  9r 


9  r 


9pfe 


(20) 
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7. 1.1.1  Computed  Observations  and  Observation  Partial  Derivatives  (con't) 


VP  ipj,  t)  =  o 


Also  one  finds  that  the  position  partial  derivatives  of  the  observation  at 
time  t  reduces  to 


VP(py,  t) 


Since  p  is  a  vector  in  inertial  space,  its  magnitude  can  be  written  as 

the  square  root  of  the  sum  of  squares  if  its  three  inertial  components  X-,  i.e. 


(P  -  P)1 


•It  ■:! 


Then  the  kzn  component  of  equation  (23)  becomes 


[vP(Py,  t)j 


It  vl‘ 


tvp 


■*  i  * 

P 


k  =  1,  2,  3). 


Using  this  result  and  that  of  equation  (22)  allows  one  to  write  equation 
(21)  as 
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7. 1.1.1  Computed  Observations  and  Observation  Partial  Derivatives  (con't 


— — —  =  •— —  .  v  (t)  [ph  =  orbit  element  or  V). 

l?|  Pfc  k 


Consider  now  the  case  when  p ^  is  a  component  of  the  observing  station's 
earth-fixed  position  vector  Then  equation  (17)  reduces  to 


or 


Since 


'  30(py,  t)‘ 

3|P| 

*E 

L  Ls 

9(?E  ) 

k  \  s)k 

3 V(pj,  t)  ' 

3 

8%  J 

k  >{%)k 

■y 

pF 

(ABCD)p 

(P  ‘  P)h 


Ap 


1  + 


GR 


lf(te}  "  W  +  A^A(te)| 


ABCD  |  r(te)-rs(tR)  +  ArA(tg) 


or 


1  + 


ApGR 


+  ArA(te))-  ^(tR) 


[  ABCD  (r(t  )  +  Ar„ (t  ) )  -  rr  (tn)  1 
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7. 1.1.1  Computed  Observations  and  Observation  Partial  Derivatives  (con't) 

where  the  ABCD  are  those  matrices  discussed  in  SECTION  2  which  transforms 
an  inertial  vector  to  the  earth -fixed  system  of  coordinates,  then  equation 
(28)  may  be  written  as 


3P(p;  t) 


(pe  •  pe)  ^  » 


3 P(P;,  t) 


EsJ  k 


|ABCD(f(te)  +  AFA(te))  -  r  (tp) 


ABCD  (r(te)  +  ArA(te))  -  r£  (tR) 

-  ^  . 

“  a* 


ABCD  (r(tg)  +  ArA(te))-  r£  (tR) 


Application  of  the  chain  rule  to  equation  (32)  gives 


3 P(P; »  t) 


\  Mp;,  t)l 


U 
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7. 1.1.1  Computed  Observations  and  Observation  Partial  Derivatives  (con't) 


Expressions  for  the  observation  partial  derivatives  with  respect  to  the  bias 
parameters  are  easily  obtained.  When  =  f^,  f^,  or  CR,  one  readily  finds 
from  equations  (17)  the  following  results: 


31 Hp/,  t) 


~r~  (t  -  tj 

'  r\ 


W(pp  t) 


J_  c 

2  f  n 


(t  -  tj) 


3 Vlpj,  t) 


7. 1.1.2  Pass  Normal  Equations  and  Parameter  Improvement 
The  results  of  the  previous  subsections  are  used  to  form  the  pass  normal 
equations.  Initially  this  is  done  using  all  of  the  observations  in  a  pass 
of  sufficient  elevation  to  form  the  residuals  (0^_C).  As  mentioned  before, 
on  subsequent  iterations  those  observations  which  have  been  identified  to  be 
inconsistent  with  the  other  observations  of  the  Dass  are  not  used  in  the 
formation  of  the  normal  equations.  An  overview  of  the  point  editor  logic 
is  presented  in  the  flowchart  of  Figure  7-1. 

The  6  and  E  matrices  are  formed  using  equations  (14),  (15),  (19),  (26),  (34), 
(35),  (36),  and  (37)  where  the  weighting  matrix  is  given  by 


Here  <j-‘  is  the  variance  for  the  (  observation;  j  is  the  Kronecker  delta; 
and  *  is  a  multiplier  which  is  set  equal  to  unity  on  the  first  iteration  and 
recomputed  for  each  subsequent  iteration  (see  SECTION  8).  For  computational 
convenience  the  B  and  £  matrices  are  partitioned  in  the  following  manner: 
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COMPUTE  ZENITH 
ANGLES  AT  OBSER¬ 
VATION  TIMES  &  TEST 
AGAINST  TOLERANCES 


MEETS 

TOLERANCE 


TAG 

OBSERVATION 


INTERPOLATE 

EOR+pkat  I 
OBSERVATION  TIMES 
(SEE  SEC  9) 


FORMyMATRIX 
EOS  1191.(261. 
(341.  (36).  (361,  (371 


FORM  OBSERVATION 
RESIDUALS  <  o-c) 
EOS  1161.  1171.  IIS) 


FORM  PASS 
NORMAL  EQUATIONS 
EOS  (131.041.(151 
1311.  139).  1401 


transform  epoch 
TO  tc<  6  TO  tcA 
REFERENCE  FRAME 
I  EOS  150)  -  155) 


SOLVE  FOR 


TRANSFORM  Ap 
TO  EPOCH  t«  & 
INERTIAL  FRAME 
EQ  1561 


APPLY 

CONVERGENCE 
TESTS 
(SEE  SEC  t! 


CORRECT 
RESIDUALS 
EO  (571 


ASSESS 
RESIDUAL 
QUALITY 
(SEE  SEC  81 


ACCEPTABLE 


TAG 

OBSERVATION 


TEST  NUMBER  OF 
REMAINING  POINTS 
i  (SEE  SEC  81 


RECOMPUTE  A 
MULTIPLIER 
EQ  (381 

(ALSO  SEE  SEC  » 


FIGURE  7-1.  Point  Editor  Process  Flow 
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Normal 

Equations 

and 

Parameter  Improvement 

(con't) 

8 

8 

8 

b  „\ 

~  e,e 

e,D 

~  e,s 

-  l 

5  D,e 

8 

D,D 

5  D,s 

-6D,b  1 

8 

8 

B 

B  . 

~  s,e 

s,D 

-  s,s 

~  s,b  | 

j 

~  b,e 

8 

b,D 

B  . 

~  b,s 

5  b,b  / 

(39) 


and 


where  the  subscripts  "e" ,  "D",  "s",  and  "b"  mean  that  the  partitioned  B 

and  E  matrices  are  formulated  using  theamatrix  of  equation  (19)  with  pfc  ; 

being  orbit  elements,  drag  parameter,  station  coordinates,  and  bias  parameters, 
respecti  vely. 

To  perform  editing  on  the  observation  residuals,  the  point  editor  first  • 

suppresses  the  submatrices  of  the  8  and  E  matrices  which  are  related  to  the 

drag  parameter  and  station  coordinates.  The  remaining  orbit  element  sub- 

matrices  are  then  re-epoched  from  the  initial  trajectory  epoch,  t  ,  to  the 

time  of  closest  approach  for  the  pass,  t^  (details  concerning  the  computation  1 

of  t are  given  in  SECTION  9).  Assuming  that  the  a  matrix  is  epoched  at 

t (i.e.  the  partial  derivatives  of  r  and  r  of  the  T  and  matrices  are 

taken  with  respect  to  the  orbit  parameters  at  t^),  one  may  use  equation  (65) 

of  SECTION  6  to  write  for  the  observation  partials  with  respect  to  the  “  .  ; 

orbit  elements: 
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7. 1.1.2  Pass  Normal  Equations  and  Parameter  Improvement  (con't) 


(^ca)* 


row 


=  (^Tc(t.)  ^Tc(t.))  *  (t.) 


(fc(t^)  fc(t.l)  J(t^)  j(tCA)  J  (tCA) 


(  50, 


th  row  ~(tCA)  2  tr. 


where 


and  V 


Using  this  result  in  equation  (14)  gives  the  desired  transformation  from  t 


~‘‘CA’=(jt0  r'(t“>  VW]T  "  [\  !_1(tcn>  5tCA'tcA> 


5(tCA) 


^ tCA^  ?trA^CA^  8  ^  ~  ^CA'  ?  tri tCA^ 

LA  ~  66  LA 


For  the  bias  partials  (i.e.  when  Pfc  =  fb,  fb,  or  CR)  one  finds  that 
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7.1. 1.2  Pass  Normal  Equations  and  Parameter  Improvement  (con't) 
so  that 


and  that 


5eb*W 


5  V 


and 


5be^W  ?be(to^  ?  ?tc/^W 


Similarly  by  applying  equation  (41)  to  equation  (15),  one  finds 

yw  =  [f'iw  5tCA^cA)  ] T  se(t0) 

and 

fb(tCA)  =  Eb(to>  ’ 


(45) 


(46) 


(47) 


(48) 


(49) 


These  results  can  be  used  to  write  for  the  "drag  parameter/station  coordi¬ 
nates  suppressed"  matrices  epoched  at  t^ 


~  (W 


(50) 

(con't  next 
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7. 1.1.2  Pass  Normal  Equations  and  Parameter  Improvement  (con't) 


?'1(tCA>  !t  'wT  §<‘o> 

LM  0c 


i  (tCA>  JtCA(tCA) 


e12  =  (tCA>  ?tCA(tCA>  T  ?eb(to> 


?be(to)  [  i  ‘ 1  tCA !  *tCA(tCA)] 


e22  =  5bb(to> 


(50) 

(con't) 


E  (tCA* 


rl<lCA>  !tCA(tCA>  I  Ee'V 


Eb'V 
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7. 1.1. 2  Pass  Normal  Equations  and  Parameter  Improvement  (con't) 


Prior  to  performing  a  solution  for  the  orbit  and  bias  parameter  corrections 
the  8(t and  E(t^)  matrices  are  transferred  to  the  reference  frame 
using  the  T  transformation  defined  in  SECTION  2.  The  submatrix  transforma¬ 
tions  are  given  by: 


-?ee(tCA) 


I  iWWl 


I  ?eb(W 


§be(W 


I 


£bb(W 


?bb(  W 


Se<W  =  I  ~Ee(  W 


fb^CA^  =  ~b^CA^ 


so  that  one  may  write 


§ee^W 


WW 


§bb(tCA)  j 


TR  82-391 


.1.1.2  Pass  Nonna!  Equations  and  Parameter  Improvement  (con't) 


Se(tCA) 


Eb^CA; 


(54) 


nd  the  associated  pass  normal  equation 


in  the  t^ 


reference  frame: 


l  I 

n  the  last  equation  AetCA  and  Ab represent  the  orbit  element  and 
ias  parameter  corrections  at  epoch  tCA  in  the  tCA  reference  frame. 


l  Cholesky  decomposition  (see  SECTION  9)  is  used  to  solve  equation  (55) 
or  the  parameter  corrections.  These  corrections  are  transformed  back  to 
he  initial  frame  at  epoch  t0  via  the  expression 


recessing  for  the  pass  is  tested  against  the  convergence  criteria  of 
ECTION  8,  and  if  they  are  met  the  updated  pass  matrices  given  by  equa- 
ions  (39)  -  (40),  along  with  both  tagged  and  untagged  observations,  are 
tored  for  processing  by  the  cross  pass  editor.  If  processing  has  not 
onverged,  the  observation  residuals  are  corrected  using  the  linear  relation 
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7. 1.1.2  Pass  Normal  Equations  and  Parameter  Improvement  (con't) 


corrected 


(0  -  C) 


current 


where  zeroes  have  been  inserted  for  drag  parameter  and  station  coordinate 
corrections,  since  none  are  computed  in  equation  (55).  The  quality  of  the 
corrected  observation  residuals  is  then  assessed  and  outlier  observations 
are  tagged  as  such.  The  differential  correction  process  then  repeats  itself 
using  only  untagged  observations  and  associated  residual s  in  equations  (38) 

-  (57). 

7.1.2  APPLICATION  TO  THE  CROSS  PASS  EDITING  PROCESS 
The  cross  pass  editor  utilizes  the  pass  matrices  (B,  E)  generated  by  the 
point  editor  to  determine  the  quality  and  useability  of  each  pass  of  obser¬ 
vations  (for  the  arc  of  interest,  i.e.  editing  span  )  in  post-PULSAR 
applications.  Those  passes  of  observations  which  do  not  meet  the  quality 
criteria  imposed  by  the  cross  pass  editor  are  flagged  as  such  during  the 
editing  process. 


The  basic  processing  of  the  cross  pass  editor  involves  first  combining  those 
pass  matrices  of  sufficient  quality  in  a  particular  manner  to  form  the  arc 
matrices  B^pp  and  E^p^  which  are  comprised  only  of  elements  related  to  the 
dynamic  parameters  (i.e.  orbit  elements  and  drag  parameter).  An  arc  solt'Mon 
for  dynamic  parameter  corrections,  Ap^p^,  is  obtained  and  used  to  compute  for 
each  pass  an  estimate  of  the  offsets  in  station  location  in  the  tr/\  reference 
frame  which  minimizes  the  observation  residuals  of  the  pass  while  keeping 
the  orbit  fixed  (hereafter,  this  will  be  referred  to  as  a  navigation  solution). 
The  quality  of  all  navigation  solutions  is  assessed  and  those  of  inadequate 
quality  are  tagged.  This  procedure  is  continued  until  the  process  converges 


with  and  E^pp  being  reformed  each  time  using  only  the  pass  matrices 
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8.1.1  ZENITH  ANGLE  EDITING  (con't) 

quality.  Furthermore,  if  the  vehicle  location  at  t^  (see  SECTION  9)  for  a 
given  pass  is  too  near  the  observing  station's  horizon,  then  it  is  assumed 
that  all  observations  for  that  pass  are  too  near  the  horizon  and  the  entire 
pass  is  identified  as  being  unusable  by  the  cross  pass  editor  and  for  ay 
post-PULSAR  processing. 

In  order  to  derive  concise  mathematical  expressions  for  the  zenith  angle  tests, 
consider  the  geometry  of  Figure  8-1  where 

N  =  Outward  directed  vector  normal  to  the  reference  ellipsoid  at 
station  location  in  inertial  reference  system, 

p  =  position  vector  of  satellite  relative  to  the  station  in  inertial 
reference  system,  and 

z  =  zenith  angle. 

One  readily  sees  that 

N 


so  that 

sin  z  =  (1  -  cos2z)^  =  [1  -  (N  •  p)2  ]2  .  (2) 

Since 

sin  z 

tan  z  =  -  ,  (3) 

cos  z 


- 1 A  i 


SECTION  8 


DATA  QUALITY  ASSESSMENT  METHODS  AND  ITERATIVE  CONVERGENCE  CRITERIA 

:  underlying  mathematics  associated  with  the  point  and  cross  pass  data  edi- 
•s  are  developed  in  detail  in  SECTION  7.  However,  little  attention  is 
'Oted  there  to  a  discussion  of  the  methods  used  to  identify  those  observa- 
ns  which  are  of  unacceptable  quality  for  use  in  post-PULSAR  processing.  Also 
itted  from  SECTION  7  is  any  discussion  of  the  criteria  used  in  PULSAR  to 
:ermine  when  the  point  and  cross  pass  processors  have  adequately  edited  the 
servation  data  and  should  terminate  their  processing,  i.e.  the  iterative 
ivergence  criteria. 

e  observation  quality  assessment  methods  and  the  iterative  convergence 
iteria  used  in  the  point  and  cross  pass  editors  are  outlined  in  some  detail 
the  following  subsections.  Also  included  below  is  a  discussion  of  the 
chnique  used  to  determine  whether  or  not  the  initial  state  vector  (and  suc- 
eding  day's  initial  conditions)  need  be  corrected  and  the  point  and  cross 
ss  editors  recycled  again  using  a  new  reference  trajectory  derived  from  the 
rrected  initial  state  vector. 

1  POINT  EDITOR  METHODS  AND  CRITERIA 

described  in  SECTION  7,  the  primary  function  of  the  point  editor  is  to 
entify  observations  on  a  pass  by  pass  basis  which  are  of  inadequate  quality, 
en  certain  conditions  are  satisfied  the  point  editor  can  also  identify  an 
tire  pass  of  observations  as  unusable  by  the  cross  pass  editor,  as  well  as 
post-PULSAR  processing.  The  methods  used  to  perform  these  identifications 
e  generally  straight  forward  and  are  described  below.  Also  included  are 
scriptions  of  the  simple  iterative  convergence  criteria  employed  by  the 
int  editor. 

1.1  ZENITH  ANGLE  EDITING 

nee  reliable  tropospheric  refraction  corrections  (see  SECTION  4)  cannot 
made  for  observations  which  lie  too  near  an  observing  station's  local 
rizon,  the  point  editor  tags  these  observations  as  being  of  inadequate 
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IJ 


and 


2,  3;  N  =  AR,  AT,  AN)  , 


(96) 


^-ARC^  ij 


where  the  a  are 


j  = 


'a  priori'  we i ght ing 


2, 


7; 


P.  = 

i 


var i ances  and  6  .  . 

'J 


e I ,  e2,  . . . ,  efi,  D) ,  (97) 
is  the  Kronecker  delta. 
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7.2  THE  PROPAGATION  OF  STATE  IMPROVEMENTS  (con't) 

Although  not  used  in  initiating  the  succeeding  day's  processing,  the  state 
covariance  8^^  is  also  propagated  to  the  succeeding  day's  epoch  for  opera¬ 
tional  convenience.  This  propagation  is  performed  in  the  following  manner: 


.ARC, 


*(ts)  !  (ts) 

n-1 

?ARC 

Hts)  1  (ts) 

.  t0  D(t0)  . 

•~to  ~D(t0)  . 

(93) 


It  should  be  mentioned  that  state  improvements  and  associated  covariance 
matrices  cannot  be  propagated  across  the  orbit  adjust  discontinuities 
which  will  exist  in  PULSAR  generated  trajectories  (see  SECTION  9).  When 
such  discontinui ties  occur,  only  the  improvements  and  covariance  associated 
with  the  state  vector  used  to  generate  the  post  orbit  adjust  trajectory  will 
be  propagated  forward  to  improve  initial  conditions  for  succeeding  days. 


7-3  'A  PRIORI1  ADDITIVE  WEIGHTING 

An  'a  priori'  additive  weighting  feature  is  generally  used  in  the  point  and 
cross  pass  editors.  The  details  of  this  feature  have  been  suppressed  for 
the  sake  of  brevity  in  the  preceeding  development.  They  can,  however,  be 
made  explicit  by  replacing:  (i)  8^  with  8^  +  everywhere  in  equations  (52), 
(58),  (60) - (62) ,  (64),  (65),  (86),  and  (88);  (ii)  each  of  the  bias  eliminated 
station  pass  matrices  with  itself  plus  Ws  in  each  summation  term  used  to  form 


-  5 

Pec-sj  inequations  (66),  (69)-(74),  (76),  and  (77);  (iii)  T'T8e 
- » T 1 


is 

T'  '  (6) 


T '  with 


1  .  -ts .  a . 

c  s  j  s  j  +  Ws)r  +  where  is  now  the  bias  eliminated  station 


pass  matrix,  everywhere  in  equations  (85)  and  (88);  and  (iv)  8  with  8  +  W 

.  .  /  \  ARC  -ARC  -ARC 

in  equations  (75),  (79),  and  (93)-  Here  the  W  matrices  are  the  'a  priori' 


weighting  matrices  given  by: 


1 


2,  3;  b. 


(94) 


(w  ).. 


-s  lj 


1 


2,  3;  x.  =  earth  fixed  coordinate), 


(95) 
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7.2  THE  PROPAGATION  OF  STATE  IMPROVEMENTS 

If  it  is  determined  that  the  initial  state  vector  needs  to  be  corrected 
using  the  arc  solution  Ap^Q  and  the  point  and  cross  pass  editing  procedure 
repeated  using  a  newly  integrated  reference  trajectory,  then  it  is  quite 
likely  that  the  initial  state  vectors  for  succeeding  days  that  were  derived 
from  the  same  initial  reference  trajectory  should  also  be  corrected  using 
the  final  arc  solution. 

Propagating  the  state  improvement  of  AP/\rc  to  a  new  epoch  time  can  be  handled 
quite  easily  by  using  the  solutions  to  the  variational  equations  that  were 
discussed  in  SECTION  6.  The  improved  state  vector  at  time  t  is  given  by 

X(t)  =  XQld(t)  +  AX(t)  ,  (89) 

where  Ay(t)  is  the  state  correction  propagated  to  the  time  t.  Using  equation 
(62) Section  6,  one  may  write  this  correction  as 


where  the  Ap^  corrections  are  the  dynamic  parameter  corrections  referenced 
to  epoch  t0.  Thus  if  the  succeeding  day's  epoch  is  t  =  t$,  then  the  improved 
initial  vector  for  that  day  is  given  by 


In  a  similar  fashion,  the  improved  state  vector  at  epoch  t0  is 

*(to>  ■  Sold(to>  +  )  *-?ARC  •  (92> 


( 


I 


-i 
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7. 1.2. 2  Navigation  and  Bias  Solutions  (con't) 


Afb  5  ?bb  (eo>  Eb(p0)  -  8be(p0)  Ape  -  Bb[)(p0)  ApD 


The  bias  solution  is  computed  from  this  expression  only  after  the  final  arc 
solution  iteration. 

A  great  deal  of  computational  convenience  is  introduced  by  using  equations 
(85)  -  (86),  since  they  can  be  written  as 

Ap  =  C  (Y)  -  C  (Y)A pp  -  C  (Y)Apn  ,  (y  =  S-,  b)  ,  (87) 

where  the  C^^Y^  are  constant  matrices  which  need  only  be  computed  once 
and  are  given  by: 

(Po)rl'1  iiTee  (p0) ;  c 

is  i  \ 

1  TiTB  (p0);  C 

s  -e 
e 

(pjrl’1  r\  (p0);  c 

-s  -  1  ~  ~  s  .0  ~ 

-ex.  -c 


=  lr?e  w 

S  -S  - 
-c  -c 


c(s^  =  T,T 

~i 


!b)  =  !bb<Eo)  Eb(P0> 

!b)  -  &<&>  &»<&>  !>  (88> 

j  °  !bb*Eo*  ?bO'?oi 
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7. 1.2.2  Navigation  and  Bias  Solutions  (con't) 
then  the  desired  result  is 


Afy  -  ?y>£o)  [!y(£o)  -  V'~o)  Ape  "  ~YD(,~o)  A~Pd] 


Equation  (84)  may  be  applied  to  obtain  the  navigation  solution  for  the  -c 
observing  station  by  letting  y  =  s ..  However,  the  cross  pass  editing  process 
performs  its  pass  data  quality  assessment  using  the  tangential  and  radial 
station  navigation  components,  AT^  and  AR^,  respectively,  in  the  t^  reference 
system.  Thus  the  above  expression  must  be  transformed  to  the  t^  system  using 
the  T'transformation  of  SECTION  2.  This  gives 


A  p  =  AT. 

~  Si 


[  T'V  (p0)  t'I  -1  [  t'TE  (p0)  -  (T\  (P0)  T)T,_1Ape 

L  ~  s  -s  •  ~  J  L  ~  s  s  -e 

JL  - C  -t  X. 


-  (p0)  nr  Vo] 

[l\  (PolT']"1  [  T,TEc  (p0)  -  T'TBc  (p0I  Ap6 

SXST  ST  Sfe 

-  l\  'Co1  ApD  ]  ’  1 


where  the  AN^  is  the  normal  station  navigation  component  and  the  subscript 
"e"  has  been  added  since  the  station  navigation  solution  is  obtained  using 
the  bias  eliminated  station  matrices.  It  should  be  noted  that  the  station 
navigations  are  obtained  for  each  iteration  associated  with  the  arc  differ¬ 
ential  correction  process. 


Bias  solutions  are  obtained  directly  from  equation  (84)  by  setting  y  =  b. 
The  expression  for  the  bias  solution  is  then 
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7. 1.2. 2  Navigation  and  Bias  Solutions 

Navigation  and  bias  solutions  are  obtained  in  the  cross  pass  editor  via  a 
linear  approximation  using  the  results  of  the  arc  solution.  To  find  this 
linear  approximation,  consider  the  normal  equation  given  by 


B  (p)  Ap 
~YY  1  -Y 


(80) 


Assume  that 


P  o  +  Ap 


ARC 


(81) 


where  p0  is  the  initial  value  matrix  for  the  parameters.  One  may  expand 
equation  (80)  in  a  Taylor  series  about  p0  to  obtain 


B  (p  ) 
~YY  ~0 


Ap 
~  Y 


s  E(P0»  + 


3  E, 


Bo 


^ARC 


(82) 


where  only  linear  terms  have  been  retained.  Since 


3Ey 


A~pARC 


t  3C 
a  1  VJ  — 

~  Y  ~  3 
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Ap 


j  9C 

2  y  'E0l  W 


ARC 
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Po 


T  3  C 

2  Y  'So1  2 
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4  D 
Bo 


a  (p0)  W  a  (pQ)  Ap 


~  Y  ~ 


fl^(pg)  W  ^[)(Pol  ApD 


?ye(Bo!  AJe 
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(83) 
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1 .1.2 A  Arc  Normal  Equation  and  Parameter  Improvement 
which  can  be  used  in  equations  (67)  -  (68)  to  obtain 

K 

8.  Ap„  +  8  Apr 
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-e  ~  e  ~c  n 

ee  eD 


"0 


E 

-C=l 


- 1 


B 


B 


es . 
a. 


s  s  • 

A.  A. 


AP«  -  6c 


's  .D 

A. 


Apn) 


and 


B 

~eDD 


a.Pd 


B" 

~e 


s  -s  . 

A.  A. 
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7. 1.2.1  Arc  Normal  Equation  and  Parameter  Improvement  (con't) 

After  the  bias  elimination  has  been  performed  upon  the  pass  matrices  they 
are  added  together,  taking  care  to  properly  segregate  the  resulting  station 
coordinate  submatrices.  For  <  observing  stations,  the  associated  bias 
eliminated  normal  equations  become 


s 
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~eee 

~£e0 
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'■Ges2 
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0 
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~£s  e 
< 

~£S  D 

< 

'£s  s 

<  K 

,  (66) 


where  the  subscript  "s  •”  means  the  matrix  is  related  to  the  x.  station. 

x. 

Thus  one  may  write 


5Se  ^  *  "ce0  *  I  K. 

X  =  1  X. 


8,  +  8  ^Pn  +  Y  8  Ap 

-'-no  ~ e  ~cnn  ~ u  4-.  -cn.  - 


U.  -I  S  •  ~r 

•  us  •  x.  D 

x.=  l  x. 


S.  Ap  *  8  Apn  +  6.  Ap  *  E.  ,  U  =1 ,  . . ,  k) .  (69) 

-‘■s-e  ~e  ~£s  -0  ~D  '"s  s.  ^s. 

X.  X.  X.  X.  X. 


Solving  equation  (69)  forAp  gives 

~  sx. 


Lpe  =  b]‘  (E_ 

~  sx:  ~°s  s .  ~"s- 

X.  X.  X. 


5  Apo 

~  e 
s  .e 
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(70) 
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7. 1.2.1  Arc  Normal  Equation  and  Parameter  Improvement  (con't) 

where  the  subscript  "a"  represents  all  but  the  bias  parameters  (i.e. 
fb,  CR)  which  are  represented  by  the  subscript  "b".  Thus  one  may  write 


8  Ap  +  8  ,  Apk 

~aot  J  a  ~ab  i  b 


?ba  ~^a  +  8hk 


~bb  Vb 


The  last  equation  may  be  solved  for  the  bias  correction  to  give 


5bb  <?b  *  5ba  V 


which  can  be  substituted  into  equation  (59)  to  give 


-1  -2 

(8  -  8  .  Bkk  8.  )  Ap  =  E  -  8  k  Bkk  Ek 

~aa  ~ab  ~bb  ~ba  ia  ~a  ~ab  ~bb  ~b 


The  last  equation  is  the  bias  eliminated  pass  normal  equation  which  can  be 
rewritten  as 


8  Ap 
~e  -la 
aa 


where  the  subscript  "e"  means  that  the  associated  matrices  are  the  bias 
eliminated  pass  matrices  defined  as 


B  8  k  8kk  8k 

-aa  ~ab  ~bb  ~ba 


E  B  .  8kk  E  k 

-  a  -ab  ~bb  -  b 
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FIGURE  7-2.  Cross  Pass  Editor  Process  Flow 
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7.1.2  APPLICATION  TO  THE  CROSS  PASS  EDITING  PROCESS  (con't) 

for  the  untagged  passes.  It  should  be  noted  that  navigation  solutions 
for  all  passes  are  obtained  during  each  arc  solution  iteration  and  their 
quality  assessed.  Any  pass  which  was  tagged  on  a  previous  iteration,  but 
is  of  adequate  quality  on  the  current  iteration  is  untagged  and  used  in 
the  arc  solution  on  the  next  iteration.  Solutions  for  station  bias  para- 
meters  (i.e.  f^,  f^,  and  C^)  are  computed  only  on  the  final  iteration.  Upon 
convergence  Ap/\pg  is  checked  to  see  if  the  initial  state  vector  should  be 
corrected  by  Ap/\Rg  and  the  process  repeated,  beginning  with  the  point  editor. 
If  so,  the  process  diverts  to  the  integrator  and  starts  again.  If  not,  the 
final  file  of  edited  tracking  data  is  output  and  the  process  terminates. 

An  overview  of  the  cross  pass  editor  process  flow  is  presented  in  Figure 
7-2.  It  should  be  pointed  out  that  should  an  orbit  adjust  occur  during  the 
editing  span,  the  cross  pass  editor  will  partition  the  span  into  two  segments 
and  cross  pass  edit  each  segment  independently.  However,  if  a  segment  is 
comprised  of  eight  or  fewer  untagged  passes,  it  will  not  be  processed  by  the 
cross  pass  editor. 

In  this  subsection,  the  mathematics  associated  with  arc  normal  equation  for¬ 
mation  and  solution,  navigation  solutions,  and  bias  solutions  are  discussed. 
The  methods  used  for  quality  assessment  and  convergence  determination  are 
developed  in  SECTION  8. 


7 . 1 . 2 . 1  Arc  Normal  Equation  and  Parameter  Improvement 
The  pass  matrices  contain  elements  which  are  station  dependent.  Since  the 
arc  matrices  must  be  formed  from  the  pass  matrices,  and  yet  be  independent 
of  station  related  elements,  one  must  devise  a  method  for  eliminating  the 
station  dependencies.  Consider  first  the  elimination  of  the  bias  elements 
from  the  pass  matrices.  Let  the  pass  normal  equation  be  partitioned  into 
submatrices  as 


6 


~ab  \  I 


(58) 


Rbrt  ?bb  /  \ 
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FIGURE  8-1.  Zenith  Angle  Geometry 
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8.1.1  ZENITH  ANGLE  EDITING  (con't) 
then  the  zenith  angle  is  obtained  from 


[1  -  (N  •  p)2  1* 


/  (N  •  p) 


To  obtain  N,  one  uses  the  transformation  from  station  geodetic  latitude, 
height  and  longitude  to  earth-fixed  cartesian  coordinates  (see  SECTION  2) 
to  find  the  earth- fixed  coordinates  of  the  station  on  the  reference  ellip¬ 
soid.  This  is  done  by  setting  the  station's  geodetic  height  above  the 
reference  ellipsoid  to  zero  giving 


(1  -  e2  sin2  <f>)^ 


cos4>  cosX 


cost))  sinX 


(1  -  e  )  sin <(> 


\  (1  '  £2)  S3  /  * 

where  <f>,  X,  and  care  the  station's  geodetic  latitude,  longitude  and  the 
eccentricity  of  the  reference  ellipsoid,  respectively. 

In  general  the  equation  for  the  surface  of  the  reference  ellipsoid  is 
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8.1.1  ZENITH  ANGLE  EDITING  (con't) 

where  be  is  the  semi-minor  axis  of  the  reference  ellipsoid.  Since  the 
eccentricity  of  the  reference  ellipsoid  is  given  by 


e 


2 


9 


then  one  may  write 

be2  =  ae2  (1  -  e2) 


Using  this  in  equation  (6)  gives 

(1  -e2)  (e2  +  f2)  +  g2  =  ae2  (1  -e2). 


(7) 


(8) 


(9) 


The  gradient  of  the  last  equation  evaluated  at  the  station  location  on  the 
reference  ellipsoid  is  an  outward  normal  to  the  surface  at  the  point  in  ques¬ 
tion.  Taking  the  gradient  of  the  last  expression  gives 


(1  -  e2)  es 


s 

1 


\ 


(1  -  £2 )  fs 


(1  -  e2)fc 


S 

2 


(10) 


\  9s  I 


where  n  is  the  normal  in  the  earth-fixed  reference  frame.  Applying  the 

—V 

ABCV  transformations  of  SECTION  2,  one  obtains  for  N: 


N  =  n 


(ID 
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8.1.1  ZENITH  ANGLE  EDITING  (con't) 

Equations  (4),  (10)  and  (11)  are  used  to  compute  the  zenith  angle  at  t^ 
and  for  each  observation  of  a  usable  pass.  The  following  two  tests  are 
applied  directly  to  the  computed  zenith  angles:  if 

<  Zjq|_  ,  accept  pass  of  observations  for  further  editing 
by  point  processor 

>  z tol  »  reject  all  observations  of  pass,  (12) 

and  if 

<  z  jq[_  accept  the  observation 

(13) 

>  2-j.qL  ,  tag  the  observation, 

where  z^  and  zt  are  the  zenith  angles  computed  from  the  reference  trajec¬ 
tory  at  t c/\  and  observation  time  t,  respectively. 

8.1.2  ORBIT  ADJUST  INTERVAL  EDITING 

If  tracking  data  falls  within  a  period  of  time  around  an  orbit  adjust  inter¬ 
val,  then  the  pass  (or  passes)  of  data  associated  with  that  tracking  data  will 
be  tagged  as  bad  data.  The  test  used  to  determine  whether  this  condition 
exists  is  given  as  follows: 

if 

tB  <  tL  (14) 

and 

t^  >  tj  ,  reject  the  pass.  (15) 
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8.1.2  ORBIT  ADJUST  INTERVAL  EDITING  (con't) 

In  the  above  inequalities  t^  and  t^  are  the  times  of  the  first  and  last 
observation  in  the  pass  and 


tB  =  t0Ag  "  ml^ 


t0Ar-  +  m2h 


where  t^  and  tp^  are  the  orbit  adjust  start  and  stop  times,  respectively, 

B  L 

h  is  the  integration  step  size,  and  the  m.  {l  =  1,  2)  are  integer  constants. 

't 

8.1.3  OBSERVATION  RESIDUAL  QUALITY  EDITING 

Two  different  types  of  quality  tests  are  performed  upon  the  observation 
residuals  in  the  point  editing  process.  These  are  called  the  absolute  and 
statistical  residual  quality  edits,  respectively.  Consider  first  the  abso¬ 
lute  residual  quality  edit.  This  is  a  coarse  check  performed  upon  all  obser¬ 
vation  residuals  in  each  pass  to  identify  those  residuals  which  are  associated 
with  observations  that  are  obviously  of  undesirable  quality.  This  test  is 
invoked  upon  the  first  iteration  of  the  point  editing  process  only  and  is 
stated  as  follows:  if 


(O^C)J 


’TOL  ’  ^ag  ^he  °t)Servat''on’ 


where  Pjql  is  a  pre-set  tolerance  value. 

The  statistical  residual  quality  edit  uses  the  weighted  least  squares  method 
to  fit  a  straight  line  through  the  observation  residuals  of  a  pass  from  which 
a  standard  error  can  be  calculated  and  used  in  the  test  procedure.  To  elaborate 
upon  this,  let 
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8.1.3  OBSERVATION  RESIDUAL  QUALITY  EDITING  (con't) 

y  =  at  +  b  (19) 


be  the  linear  function  of  time  t  which  one  wishes  to  fit  to  the  observations. 
Using  the  results  of  SECTION  7  one  finds  that  the  required  normal  equation 
may  be  written  as 


where  N  is  the  number  of  observation  residuals  in  the  fit;  X  is  the 
multiplier  defined  by  equation  (38)  of  SECTION  7;  and  is  the  time 
of  the  ^th  observation.  Matrix  equation  (20)  is  solved  for  a  and  b  which 
are  used  with  equation  (19)  to  compute  the  standard  error  for  the  fit  using 
the  expressions: 

r  .  =  (0  -  C)  .  -  y .  ( <  =  1,  2,  3,  . . . ,  N) ,  and  (21) 

^  -t  corrected 


(22) 


Equation  (22)  is  used  to  perform  the  statistical  residual  quality  edit 
on  each  point  editor  iteration  but  the  first.  This  test  is  stated  as 
follows:  if 


I  ri  I  >  2-5  o 


,  tag  the  observation  . 


(23) 


8.1.3  OBSERVATION  RESIDUAL  QUALITY  EDITING  (con't) 

As  noted  in  SECTION  7,  the  tagged  observations  are  used  in  each  iteration  to 
perform  the  linear  fit  and  in  equations  (21)  -  (23).  They  are  not  used  in 
the  formation  of  the  pass  normal  equations. 

The  X  multiplier  in  equation  (20)  (and  equation  (38)  of  SECTION  7)  is  set 
equal  to  unity  on  the  first  iteration  and  is  computed  on  successive  iterations 
from  the  expression 


V  =  (S/N)  X 

new  v  p  previous 


(24) 


where 


(S/N)p  * 


V-  -T 


e  e 

~  s  ~s 


/( 


MAX  (1,  Nr 


Ns) 


(25) 


In  this  equation  Nq  is  the  number  of  untagged  points  in  the  pass;  Ns  is  the 
number  of  best  determined  parameters  used  in  the  solution  to  the  pass  normal 
equations;  e <.  is  a  transformed  E  matrix  used  in  the  solution  to  the  pass 
normal  equations;  and  V  is  the  variance  given  by 


V 


NG 

E 

i=\ 


(0  -  C)  . 

±  corrected 

Xa  .2 


(26) 


The  N$  number  and  the  matrix  are  discussed  in  the  Cholesky  decomposition 
subsection  of  SECTION  9. 

8.1.4  OBSERVATION  QUANTITY  EDITING 

As  the  point  editor  iterates  during  its  processing,  it  will  tag  observations  as 
being  of  insufficient  quality.  It  may  also  untag  previously  tagged  observa¬ 
tions  if  they  return  to  the  quality  bounds  established  during  each  iteration 
(see  SECTION  8.1.2).  This  processing  may  produce  a  situation  where  too  few 
observations  remain  in  a  pass  for  a  successful  solution  to  the  associated  pass 


8.1.4  OBSERVATION  QUANTITY  EDITING  (con't) 

normal  equation  to  be  obtained.  Such  passes  are  then  rejected.  The  standard 
used  to  determine  this  is  straight-forward  and  given  as  follows:  if 


NR  <  .01  P  •  N  ,  reject  the  pass,  (27) 

where  NR  is  the  number  of  untagged  observations  remaining;  P  is  the  pre¬ 
determined  percentage  of  original  observations  which  must  be  untagged  and 
below  which  the  pass  is  rejected;  and  NQ  is  the  original  number  of  observa¬ 
tions  in  the  pass. 

8.1.5  ITERATIVE  CONVERGENCE  CRITERIA 

Observation  data  is  processed  by  the  point  editor  on  a  pass-by-pass  basis. 

When  the  last  pass  has  been  processed,  the  point  editing  procedure  is  terminated. 
This  termination  is  keyed  from  an  end  of  file  condition. 

In  order  that  the  point  editor  exit  at  the  proper  time  from  the  iterative 
loop  for  a  pass,  one  of  the  following  two  conditions  must  be  satisfied: 
if  either 


Ni 


MAX, 


or 


,  exit  iteration  loop  normally  and  (28) 
process  next  pass. 


where  Nr  is  the  number  of  iterations  performed  for  a  pass; 

P 

maximum  number  of  iterations  allowed  for  a  pass;  and  Nq^  is 
untagged  observations  for  iteration  fe. 


nIMAXp  is  the 
the  number  of 


8.2  CROSS  PASS  EDITOR  METHODS  AND  CRITERIA 

The  cross  pass  editor  uses  the  statistics  of  the  navigation  solutions  for 
all  passes  to  identify  those  passes  of  observation  data  which  lack  the 
quality  required  for  post-PULSAR  processing.  The  methods  used  to  identify 
those  passes  are  described  below.  The  iterative  convergence  criteria  used 
by  the  cross  pass  editor  are  also  given.  It  should  be  mentioned  that  if  an 
orbit  adjust  occurs  during  the  editing  span,  then  the  tests  and  criteria 
discussed  below  are  applied  independently  to  each  trajectory  segment  (refer 
to  subsections  7.1.2  and  9.3.4). 


8.2.1  NAVIGATION  SOLUTION  QUALITY  EDITING 

In  order  to  identify  those  passes  of  data  which  are  outliers,  on  each  itera¬ 
tion  the  cross  pass  editor  performs  a  coarse  tolerance  test  to  identify  and 
to  tag  those  passes  which  have  tangential (AT  ■)  or  radial  (AR.)  navigation 
solutions  (see  equation  (85)  of  Section  7)  that  are  obvious  outliers.  The 
remaining  untagged  AT^  and  AR^  are  used  to  fit  a  quadratic  equation  from  which 
a  standard  deviation  is  computed  and  used  to  edit  the  AT^  and  AR^  for  the 
remaining  untagged  passes. 

To  be  more  explicit,  the  coarse  tolerance  test  used  by  the  cross  pass  editor 

f  h 

is  as  follows:  if  for  the  V  pass 


lARJ  5  nT0L 


AR, 


or 


lATJ  >  nT0L 


ATc 


tag  the  pass. 


(29) 


where  the  Hjql  are  predetermined  tolerance  values.  The  remaining  untagged 
AR^  and  AT^-  are  fit  with  a  quadratic  of  the  form 
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8.2.1  NAVIGATION  SOLUTION  QUALITY  EDITING  (con't) 

The  required  normal  equation  may  be  written  as 


where  N  is  the  number  of  untagged  navigation  solutions  and  8^  is  either 
AR^-  or  AT^,  depending  upon  which  is  being  fit.  Using  the  definition 


r . 

/C 


AR  l 

or 


1»  2,  . . . ,  N) , 


the  test  u-ed  to  quality  edit  the  AR  ■  and  AT.  navigation  solutions  may 

'C  'C 

be  stated  as  follows:  if  for  the  jth  arc  solution 


r.2 


,  tag  the  pass. 
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8.2.1  NAVIGATION  SOLUTION  QUALITY  EDITING  (con't) 
In  this  expression 


N 


1=  1 


2  X.  L_ 

and  a;  ^ is  the  standard  error  for  the  (j  -  l)in  arc  solution  given  by 


°M 


N 

I 

i-l 

10s 


(%  -  yj)2: 


j- 1 


for  j  >  1 

for  j  =  1. 


(35) 


In  the  last  two  equations  ^  represents  AR^  or  AT^,  depending  upon  which  is 
being  edited,  and  N  is  the  number  of  untagged  passes.  This  process  (i.e. 
equations  (30)  -  (35))  is  repeated  m  times,  where  m  is  a  integer  constant. 
Each  of  the  m  times  this  process  is  repeated  only  untagged  passes  are  used 
to  form  the  test  conditions  (i.e.  (30),  (31),  (34),  (35)),  but  both  tagged 
and  untagged  passes  are  tested  using  inequality  (33).  Those  passes  which  are 
tagged  after  the  m^*1  iteration  are  omitted  from  the  arc  solution.  However, 
during  each  arc  solution  iteration,  those  passes  which  were  tagged  are  used 
again  to  generate  new  navigation  solutions  and  their  quality  again  assessed. 
If  it  is  adequate  the  tags  are  removed  and  the  passes  are  used  in  generating 
succeeding  arc  solutions. 


8.2.2  ITERATIVE  CONVERGENCE  CRITERIA 

In  order  that  the  cross  pass  editor  exit  from  its  iterative  arc  solution  loop 
at  the  proper  time,  one  of  the  following  two  conditions  must  be  satisfied: 
if  either 


MAX. 


or 


5  . 
1 


exit  iteration  loop  normally, 


(36) 
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8.2.2  ITERATIVE  CONVERGENCE  CRITERIA  (con't) 

where  Nj^  is  the  number  of  iteration  loops  completed  by  the  cross  pass 
editor;  nImaxp^s  t*ie  max’mum  number  of  iterations  allowed  for  the  cross 

v 

pass  editor;  and  is  the  set  of  passes  that  were  tagged  on  iteration 

k. 

8.3  STATE  VECTOR  VALIDITY 

As  an  option,  the  arc  solution  A p/\pc  obtained  from  the  last  iteration  of 
the  cross  pass  editor  can  be  tested  to  see  if  the  initial  state  vector 
should  be  corrected  to  produce  an  improved  reference  trajectory.  If  the 
testing  proves  positive  then  the  state  vector  at  trajectory  epoch  t0  is  cor¬ 
rected  using  the  results  of  SECTION  7.2;  a  new  reference  trajectory  is  gen¬ 
erated;  and  the  point  and  cross  pass  editing  processes  are  repeated.  It 
should  be  mentioned  that  when  an  orbit  adjust  occurs  during  the  editing  span, 
this  process  is  repeated  independently  for  each  trajectory  segment. 

In  order  to  devise  a  method  for  determining  the  acceptability  of  the  initial 
state  vector  consider  the  pass  variance  defined  by: 

N 

=  (q_^c(p0))T  u  «U_c(p0))  ,  (37) 

' 

x>l 

where  N  is  the  number  of  observations  during  the  pass  and  pQ  is  the  set  of 
parameters  to  be  corrected.  Suppose  the  arc  solution  Ap^^  has  been  obtained 
upon  convergence  of  the  cross  pass  editor.  If  this  arc  solution  is  used  to 
correct  the  initial  state  vector  and  regenerate  a  new  reference  trajectory, 
then  some  new  A p  will  be  introduced  into  the  pass  normal  equations.  The 
associated  pass  variance  then  becomes  V(p0  +  A p)  which,  when  expanded  in  a 
Taylor  series  about  p0,  has  the  form 


-153- 


8.3  STATE  VECTOR  VALIDITY  (con’t) 


V(p0  +  Ap)  =  V(P0)  + 


T 

po 


(38) 


Since 


then 


-  2E 


3E 

-  2  —  =28 
9  P 


so  that  one  may  write 


(39) 


(40) 


V(P0  +  Ap)  ~  V(P0)  "  2E(pQ)TAp  +  ApT  B(p0)  Ap  . 


(41) 


Using  equation  (58)  of  SECTION  7,  the  last  equation  may  be  rewritten  as: 


V(Po  +  Ap)  =  V(p0)  -  2^(P0^APa  '  2!h(po)Apb  +  Ap  A~Po 


+  ApI  ?(po)Apb  +  AA  ?(po)Apa  +  Apk  5(po)Apb  ’  (42) 


a  ~  . ~  V  ~  b  ~u’  ~  a  ~  b  i» 

ab  ba  bb 


where  as  before  subscript  "a"  represents  all  dynamic  and  station  location 
parameters  and  subscript  "b"  represents  the  bias  parameters.  If  equation 
(61)  of  SECTION  7  is  substituted  for  Ap^  in  the  last  equation,  one  finds 
that  after  some  algebraic  manipulation: 


8.3  STATE  VECTOR  VALIDITY  (con ' t) 


V(Po  +  Af)  “  V(P0)  +  APa  |8aa(p0)  -  ^b(p0)  Bbb(p0)  ?ba(P0)JAf 


2  I  ( Pn )  -  E^PJ  8  !(P„)  8  (Prt)  I  AP„ 


-a 


*b  ~°  ~bb  ~°  ~ba 


-o  ~  a 


£<&>  Eb<e0> 


By  inspection  it  is  seen  that  the  two  terms  in  brackets  in  the  last  equa¬ 
tion  are  the  bias  eliminated  pass  matrices  given  by  equations  (64)  -  (65) 
in  SECTION  7.  Thus 


V(p0  +  A p)  «  V(p0)  +  ApT  B£  (pQ)  Apa  -  2  E£  (p0)  Ap 

~  a  aa~  a 


The  last  equation  can  be  summed  over  all  passes,  taking  care  to  properly 
segregate  the  resulting  station  coordinate  submatrices.  This  gives 


V'(P0  +  Ap) 


V'(p0)  +  Ap’T  8'(p0)  Ap'  -  2  E'T(p0  )Ap '  - 


-  l 


all 

passes 


b  (p0: 

’bb  ~ 


(P0} 


where  Ap',  8'(p0),  and  E'(p  )  have  the  form  given  by  matrix  equation  (66) 
of  SECTION  7.  Station  elimination  can  be  performed  using  equation  (70)  of 
SECTION  7  to  yield  the  final  result: 
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TABLE  9-1.  TABULATION  OF  COWELL  COEFFICIENTS,  C -  ,  AND  ADAMS  COEFFICIENTS  a- 


0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 


C  J 

1/12 

0 

-1/240 

-1/240 

-221/60480 

-19/6048 

-9829/3628800 

-407/172800 

-330157/159667200 

-24377/13305600 

-4281164477/2615348736000 

-70074463/47551795200 

-1197622087/896690995200 

-97997951/80472268800 


1/2 

-1/12 

-1/24 

-19/720 

-3/160 

-863/60480 

-275/24192 

-33953/3628800 

-8183/1036800 

-3250433/479001600 

-4671/788480 

-13695779093/2615348736000 

-2224234463/475517952000 

-132282840127/31384184832000 
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1.3.1  RUNNING  PROCEDURE  (con't) 


rz(h2in)  =  v"2(h2tn_1)  +  rx(h>'iu)  (27) 

and 

V-1  (At2  xn)  =  +  V°(^2?n)  ,  (28) 

where  V~2hzxn  and  V_1fi2xn  are  called  the  second  and  first  sums,  respec¬ 
tively.  They  are  computed  during  the  starting  procedure  and  are  extrapolated 
forward  as  required  during  the  running  process. 


These  extrapolated  differences  are  used  to  predict  the  position  and  velocity 
at  t  i  through  the  application  of  the  predictor  formulae 


n+1 


=  7  (/i2 


V 


I  C; 


Vl> 


/=0 


(29) 


and 


n+1 


1 

h 


2 


I 

i=o 


7J(fi 


2  x  ) 
xn+l' 


(30) 


where  m  is  the  order  of  integration,  i.e.  the  number  of  acceleration  values 
used  in  the  predictors,  and  the  C j  and  a'-  are  the  Cowell  and  Adams  coeffi¬ 
cients,  respectively.  These  coefficients  are  tabulated  in  Table  9-1.  These 
predicted  position  and  velocity  values  are  used  in  Equation  (1)  of  SECTION  5 
and  equation  (7)  of  SECTION  6  to  produce  the  predicted  accelerations 

h2  x„, ,  .  An  extrapolated  acceleration  h2x  ,,  is  obtained  from 

n+1pred.  n+1extrap. 

equation  (26)  when  k  =  0  and  the  two  are  differenced  to  obtain  the  correction 


6(h2x)  =  h: 


n+1 


pred. 


extrap. 


(31) 


-168- 


TR  82-391 


9.3  THE  COWELL  INTEGRATOR  (con't) 

trajectory  epoch  so  that  the  eight  point  interpolation  method  described 
earlier  can  be  used  in  the  neighborhood  of  the  start  time.  The  following 
subsections  describe  the  basic  algorithm,  the  starting  procedure,  and  those 
special  procedures  mentioned  above. 

9.3.1  RUNNING  PROCEDURE 

Once  the  starting  procedure  has  been  completed  (see  next  subsection)  the 
integrator  operates  in  the  running  mode.  The  integration  step  from  tn  to 
tn+1, where 

tn+l  =  tn  +  k  ,  (22) 

and  h  is  a  fixed  integration  step  size,  is  initiated  by  extrapolating  the 
backward  difference  table  forward  one  time  step.  The  backward  differences 
V^k2x ( t^ )  are  defined  by  the  following  relationships: 

V(fi2xn)  =  v(fr2x(tn))  =  h2[x(tn)-  x(tnl)J  (23) 


Vk(h2xnJ  =  V  [  Vk'l(k2  xj  ]  .  (24) 

(Here  x  will  be  used  to  represent  both  r  and  Yp^,  since  both  the  equations 
of  motion  and  the  variational  equations  are  integrated  simultaneously  using 
the  same  algorithm).  The  extrapolation  is  performed  as  follows  by  holding 
the  highest  difference  fixed: 


Ah2  xn) 

(25) 

■ 

Vfe(h2xn)  +  Vfe+1U2  xn+1). 

k  =  m  -  1 ,  • 

o 

(26) 

4 


- .  -  j 
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9.2.2  COMPUTATION  OF  ECCENTRIC  ANOMALY  (con't) 


This  result  can  then  be  used  with  equation  (16)  to  construct  the  g(E') 
function  and  the  following  recursion  relation: 


E1  =  E  1 
L  n+1  n 


En’  -  e"  sinEn‘  -  C 


1  -  e"  cosEn' 


The  eccentric  anomaly  is  computed  by  repeating  the  calculation  of  the  last 
equation  until 

|E*  ,  -  E  *1  <  1CT8  radian  .  ( 

1  n+i  n'  — 


When  this  condition  is  satisfied,  then 


E' 


E 


n+1 


( 


The  computations  of  equation  (18)  are  initiated  using 


E 


1 


l"  + 


( 


e"  sinT  \ 
1  -  e"  cosf"  / 


9.3  THE  COWELL  INTEGRATOR 

The  Cowell  predictor-corrector  method  is  used  by  PULSAR  to  simultaneously 
numerically  integrate  the  satellite  equation  of  motion  (see  equation  (1)  of 
SECTION  5)  and  the  associated  variational  equations  (see  equation  (7)  of 
SECTION  6).  This  numerical  integrator  is  a  fixed  step/mul tisteD  algorithm 
based  upon  the  integration  of  a  backward  difference  formula  fitted  to  a  set 
of  acceleration  data.  It  is  not  a  self  starting  process  and  uses  a  unique 
technique  for  supplying  a  backward  difference  table  at  the  trajectory  start 
time  necessary  for  operating  in  the  running  mode.  Furthermore,  special  pro¬ 
cedures  are  required  for  trajectory  generation  in  the  neighborhood  of  orbit 
adjust  intervals,  as  well  as  for  generating  trajectory  data  prior  to  the 
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9.2.1  COMPUTATION  OF  TIME  OF  CLOSEST  APPROACH  (tCA)  (Con’ t) 

or  until  the  maximum  iteration  count  is  reached.  If  equation  (13)  is  satis¬ 
fied  then 


However  if  the  maximum  iteration  count  is  reached,  then 


est 


where  test  is  the  estimated  t ^  supplied  for  the  pass  by  the  observing 
station.  This  is  obtained  by  sensing  a  change  of  sign  in  the  range  differ¬ 
ence  measurements.  Should  no  sign  change  occur,  t  ^  is  set  by  the 
station  to  the  time  of  the  second  observation.  It  is  possible  that  this 
will  be  a  frequent  occurrence  for  the  satellites  supported  by  the  PULSAR 
system,  since  their  continuously  varying  attitudes  can  result  in  partial 
passes  of  observation  data  that  may  contain  no  t^^  which  satisfies  equation 
(10).  It  should  be  mentioned  that  equation  (12)  for  i  =  1  is  initiated 
using  tj  -  te$t. 

9.2.2  COMPUTATION  OF  ECCENTRIC  ANOMALY 

As  mentioned  in  SECTION  5,  the  Newton  -  Raphson  iteration  is  used  to  compute 
the  eccentric  anomaly  associated  with  the  mean  anomaly  and  eccentricity  gen¬ 
erated  by  the  Brouwer  orbit  predictor.  This  is  done  by  iteratively  solving 
the  Kepler  equation  given  by  equation  (72)  of  SECTION  5.  One  rewrites  this 
equation  as 

f(E’)  =  E*  -  e"  sinE 1  -  f"  =  0  ,  ( 

where  E 1 ,  #■"  and  e"  are  the  eccentric  anomaly,  the  mean  anomaly,  and  the 
eccentricity,  respeci tvely .  This  expression  may  be  differentiated  with 
respect  to  E'  to  give 


f'(E')  = 


d  f(E') 


1  -  e"  cosE' 


9.2.1  COMPUTATION  OF  THE  TIME  OF  CLOSEST  APPROACH  (tcA) 


where 

•  • 

P  =  r  - 

p  =  r  - 
R  =  a  x 
and 

R  =  w  x 

The  satellite  velocity  and  acceleration  vectors  are  interpolated  using 
equations  (2)  and  (3).  The  angular  velocity  of  the  earth  is  given  by 
equation  (18)  of  Section  5. 

Equations  (11)  are  used  to  construct  the  g(t)  function  and  the  associated 
recursion  relation: 


R 

-»- 

R 

*> 

R 

R 


t-,,  =  t- 

•C+l  ^ 


p(t ; )  •  p(t,) 


p(t^)  •  p(t^)  +  p(t^)  •  p(t^) 


To  compute  t^  this  recursion  relation  is  repeated  until 
|tf+i  '  <  *05  second 
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OBSERVING 

^STATION 


^SATELLITE 


EARTH 


r|t)  =  INERTIAL  POSITION  VECTOR  OF  SATELLITE  AT  TIME  t 

3(t)  =  INERTIAL  POSITION  VECTOR  OF  STATION  AT  TIME  t 

pit)  =  SATELLITE  POSITION  VECTOR  RELATIVE  TO  THE  STATION  AT  TIME  t 

pit)  =  SATELLITE  VELOCITY  VECTOR  RELATIVE  TO  THE  STATION  AT  TIME  t 


FIGURE  9-1.  Satellite-Observing  Station  Relative  Geometry 


TR  82-391 


9.2  NEWTON  -  RAPHSON  ITERATION 

The  Newton-Raphson  method  is  an  iterative  technique  which  is  used  to  find  the 
zeroes  of  adifferentiable  function.  It  is  used  in  PULSAR  to  find  times  of 
closest  approach  (t^)  for  individual  passes,  as  well  as  to  solve  Kepler's 
equation  for  the  eccentric  anomalies  of  NAVSAT  satellites  at  a  specific  time 
for  the  purpose  of  clock  calibration  (see  SECTION  3).  In  general  the  method 
may  be  stated  as  follows:  if  f(x)  is  a  function  such  that  f'(x)  = 


d  f(x) 


t  0,  and  pQ  is  an  approximation  to  a  solution  of  f(x)  =  0,  then  the 


sequence  {pw}  which  converges  to  the  zero  at  x  =  p  can  be  generated  recursively 
by  the  relation  pn  =  g(p  ^),  n  =  1,  2,  ••*,  where 


f  ( x ) 
f'(x) 


The  application  of  this  method  to  finding  t^  and  solving  the  Kepler  equa¬ 
tion  are  discussed  below. 

9.2.1  COMPUTATION  OF  TIME  OF  CLOSEST  APPROACH  (^CA) 

The  time  of  closest  approach  of  a  satellite  to  an  observing  station  is  the 
point  in  time  at  which  the  relative  velocity  component  along  the  line  of  sight 
changes  in  sign.  In  order  to  formul ate ^ this  in  a  precise  way  consider  the 
geometry  of  Figure  9-1,  where  p(t)  and  p(t)  are  the  satellite  position  and 
velocity  vectors  relative  to  the  observing  station,  respectively.  By  inspec¬ 
tion  one  sees  that  the  t^  is  the  solution  of  the  equation 

p(t)  .  p(t)  =  0  .  (1C 

To  apply  the  Newton-Raphson  method  to  solving  this  equation  for  t^,  one 
takes  the  derivative  with  respect  to  time  of  the  last  equation  giving 


p(t)  •  p(t) 


>(t)  •  p(t)  +  p(t)  •  p(t)  , 
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9.1  LAGRANGE  INTERPOLATION 

Numerical  solutions  to  the  variational  equations,  satellite  ephemerides,  sun 
and  moon  ephemerides,  nutation  angles,  pole  wander  values,  and  UT1  -  UTC 
time  corrections  are  stored  at  discreet  equally  spaced  times  for  use  within 
PULSAR.  In  order  to  obtain  these  data  at  times  other  than  those  stored,  it 
is  necessary  to  employ  an  interpolation  procedure.  PULSAR  utilizes  the 
Lagrange  interpolation  method  of  various  orders  to  obtain  these  data.  This 
method  involves  fitting  polynomials  to  the  data  points  being  interpolated 
and  centering  the  fits  such  that  equal  numbers  of  points  lie  on  either  side 
of  the  interpolation  time. 

The  polynomials  for  interpolating  each  type  of  data  stored  by  PULSAR  are 
of  the  form: 

n 

»<t)  *  I  X  <‘>»  <V  • 

-t=0 

First  time  derivatives  of  satellite  position  and  variational  equation  solu¬ 
tions  and  second  time  derivatives  of  satellite  position  are  computed  using: 
n 

- 1 X 


i/(t)  =  Y.  My  (*t)  • 

1= 0  * 

In  the  above  expressions  t  is  the  interpolation  time,  t.  are  the  times  for 
which  the  data  to  be  interpolated  exist,  and 
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SECTION  9 
NUMERICAL  METHODS 

In  the  previous  sections  a  mathematical  description  of  the  PULSAR  system 
has  been  developed  and  discussed.  However,  since  the  system  is  largely 
automated,  a  gap  still  exists  between  the  purely  analytical  treatment  of 
the  preceeding  text  and  the  practical  implementation  of  the  analytics  into 
the  framework  of  the  computational  capability  of  tne  electronic  computer. 

For  example,  the  analytic  development  provides  equations  of  motion  which 
must  be  numerically  integrated  by  the  computer  to  generate  a  usable 
ephemeris;  but  the  analytical  development  does  not  prescribe  a  method  by 
which  the  computer  can  do  this.  It  is  the  purpose  of  this  section  to  close 
this  gap  by  describing  in  some  detail  those  special  numerical  techniques  em¬ 
ployed  by  the  PULSAR  system  for  use  on  the  computer. 

Presented  in  the  following  subsections  are  overviews  of  the  numerical  tech¬ 
niques  used  by  PULSAR  for: 

U)  data  point  interpolation; 

(<l)  the  solution  of  equations  by  fixed-point  iteration; 

ilii)  numerically  integrating  the  equations  of  motion  and  variational 
equations; 

iiv)  solving  the  pass  normal  equations;  and 
(v)  compacting  observation  data. 

Should  further  detail  concerning  these  methods  be  required  the  reader  is  refer 

i  z 

red  to  several  excellent  numerical  analysis  texts.  ’ 

1  Blum,  E.  K. ,  Numerical  Analysis  and  Computation  Theory  and  Practice,  Addison 
Wesley  Publishing  Co.,  1972. 

2  Moursund,  0.  6.  and  Duris,  C.  S.,  Elementary  Theory  and  Application  of 
Numerical  Analysis,  McGraw-Hill  Book  Company,  1967. 
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8.3  STATE  VECTOR  VALIDITY  (con't) 

Jj  In  the  above  inequalities  is  the  number  of  times  the  editing  cycle 

has  been  repeated;  the  maximum  number  of  such  cycles  allowed; 

and  6^  U  =  1,  2)  are  predetermined  tolerance  values. 

|  As  noted  in  SECTION  7.2,  if  a  recycling  is  required,  the  state  vector  for 

succeeding  days  should  be  improved  using  the  correction  propagation  techni¬ 
ques  discussed  in  that  section. 

r — 

Ml 


TR  82-391 


8.3  STATE  VECTOR  VALIDITY  (con't) 

where  (S/N)*.^  is  the  station  and  bias  reduced  arc  signal  to  noise  given 
by: 


(S/N), 


V(B0)  "Z!~Gr  ^Bo^-e  ^Bo^  ~bb^~0^  ~b^0^ 


k  Sfe 


k  al  1 

asses 


From  this  one  may  also  compute  the  expected  percent  change  in  the  arc  variance 
due  to  a  correction  ApARC  applied  to  the  initial  state  vector.  This  is  given 
by  the  expression 


A(S/N); 


<S/N4r  -  ^’predicted 
(S/N)l„ 


These  results  are  used  to  formulate  the  criteria  used  by  PULSAR  to  determine 
if  a  new  reference  trajectory  need  be  generated  and  another  editing  cycle 
performed.  If  this  option  is  selected,  then  the  criteria  may  be  stated  as 
follows:  if 


NCYCLE  <  NCYCLE» 


(S/N) 


predicted 


>  i  i 


then  correct  the  initial  state  vec¬ 
tor,  regenerate  reference  trajectory, 
and  recycle  through  point  and  cross 
pass  editors. 
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8.3  STATE  VECTOR  VALIDITY  (con't) 

v'(p0  +  ap)  =  v'ip0>-y^  <p„> (p0> ?e  <fo> 

t  sfe  sfesfe  sfe 

■  I  ETb(fio>  !bb(Po>  !b<Eo> 
all 

passes 

Thus  an  estimate  of  the  variance  for  an  orbit  with  corrected  parameter  values 
;i  +  i,p  can  be  predicted  by  using  the  variance  for  the  uncorrected  orbit  and 
reducing  it  by  station,  bias,  and  arc  solution  terms. 

To  formulate  in  a  concise  way  the  two  criteria  used  to  determine  if  the  tra¬ 
jectory  must  be  regenerated  from  an  improved  state  vector  and  the  editing 
cycle  repeated,  one  further  definition  is  required:  define  the  arc  signal 
to  noise  ratio  (S/N)  as 


E  (p„ )  Ap 
-ARC  ~°  -  ARC 


(46) 


(S/N)2 


V' 

M 

I", 


(47) 


■/here  N(  is  the  number  of  observations  associated  with  the  l pass  and  M 
is  the  number  of  passes.  Then  equation  (46)  can  be  used  to  obtain  a  pre¬ 
diction  of  the  signal  to  noise  value  for  a  corrected  orbit,  i.e. 


;/n) 


predicted 


(S/N) 


SbR 


EARC(£o)  Aparc 


(48) 
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9.3.1  RUNNING  PROCEDURE  (con't) 

This  difference  is  then  used  to  correct  the  position  and  velocity  by  appli¬ 
cation  of  the  expressions: 


n+1 


ln+l 


corrected 


corrected 


=  x  +  5x 
n+1  n+1 


xn+l  +  5Vl  ’ 


(32) 


(33) 


where 


*  x 


n+1 


=(Z  cj) 


(34) 


J= 0 


and 


ex 


n+1 


-  -MI  ■v)6(’,2j) 


(35) 


f =0 


Equations  (32)  -  (33)  are  used  to  compute  new  accelerations  and  this  process 
is  continued  until  the  number  of  such  iterations  exceeds  that  of  an  in out 
control  parameter.  The  extrapolated  difference  table  is  corrected  using  the 
expression 


VJ(h‘ 


n+1 


corrected 


=  VJ(h 


n+1 


extrapolated 


Ah 


(36) 


for  j  -  0,  1,  •••,  m.  This  completes  the  integration  from  tn  to  t  +j.  The 
complete  process  is  repeated  for  tn+^,  tn+3>  etc.  A  diagram  is  presented  in 
Figure  9-2  which  shows  in  flow  form  the  running  procedure  (as  well  as  the 
starting  procedure  described  below). 
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9.3.2  STARTING  PROCEDURE 

The  Cowell  integration  process  is  initiated  by  using  the  initial  conditions 

•  • 

at  trajectory  epoch  xQ  =  x(tQ)  and  xQ  =  x(tQ)  to  compute  the  accelerations 
h2  xQ  from  equation  (1)  of  SECTION  5  and  equation  (7)  of  SECTION  6.  The 
difference  table  is  initialized  to  zero  for  fe  >  0,  i.e. 


Vk(h2 


0  ,  k  =  1,  2,  3,  • • • ,  m. 


(37) 


and  the  first  and  second  sums  are  computed  by  transposing  equations  (29)  - 
( 30 ) ,  i.e. 

m 

v_1(h2x_i)  =  ktQ  -  ^  a'j  Vj\h2  xQ)  (38) 

i= 0 


and 


-  2 

V 


m 


ZcyvVv  ■ 
i=o 


(39) 


Integration  is  then  performed  for  m  time  steps  as  described  in  the  previous 
subsection  for  the  running  procedure,  except  that  the  corrections  to  x  and 
x  are  given  by 


'(£cj) 

y=o 


and 


6{h7  x) 


(40) 


(41) 
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9.3.2  STARTING  PROCEDURE  (con 1 1) 

and  only  the  first  p  differences  are  corrected  by  5{h  x)  ,  i.e.  equation 
(36)  applies  for  j  =  0,  1,  ■  • ,  p< m.  After  the  mth  time  step  has  been  completed, 
the  difference  table  is  extrapolated  backward,  holding  the  mth  difference  con¬ 
stant,  to  obtain  an  improved  difference  table  at  the  trajectory  epoch  using 
equations  (25)  -  (26). 

This  completes  the  first  cycle  of  the  starting  procedure.  Cycling,  which  starts 
with  the  computation  of  the  first  and  second  sums  using  equations  (38)  -  (39), 
is  repeated  until  the  conditions 


(42) 

(43) 


are  met  for  all  p,  where  the  U  =  1,  2)  are  tolerance  values.  The  starting 
procedure  is  then  terminated  and  the  running  procedure  is  invoked.  It  should 
be  noted  that  during  the  first  m  steps  of  the  running  mode  the  procedure  for 
the  starting  mode  is  followed  in  that  equations  (40)  -  (41)  apply  and  only  the 
first  p  differences  are  corrected.  Figure  9-2  shows  the  flow  of  the  starting 
and  running  procedures. 

9.3.3  BACKWARD  INTEGRATION 

In  order  to  apply  the  eight  point  Lagrange  interpolation  method  in  the  vicinity 
of  the  trajectory  epoch  t0  it  is  necessary  to  integrate  the  satellite  equations 
of  motion  (equation  (1)/SECTI0N  5)  and  the  variational  equations  (equation  (7)/ 
SECTION  6)  for  several  time  steps  backward  in  time  from  tQ.  It  should  be  noted 
that  the  integration  should  also  extend  several  time  steps  beyond  the  trajectory 
stop  time  t  so  that  the  interpolation  technique  can  be  applied  in  the  neighbor¬ 
hood  of  te.  The  latter  case  presents  no  problem,  since  the  normal  integration 
process  described  in  the  previous  subsections  can  be  continued  for  some  At 
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9.3.3  BACKWARD  INTEGRATION  (con't) 

beyond  t  .  However  the  backward  integration  case  must  be  treated  in  a  special 
fashion. 

The  backward  integration  process  for  integrating  the  equations  of  motion  is 
initiated  by  changing  the  initial  conditions  from  (r(tQ),  r(tQ)j  to  ^ r(tQ) ,  -  r(tQ) )  . 
This  sets  the  satellite's  motion  in  the  proper  direction  as  the  time  is  stepped 
backwards.  However  in  order  to  produce  the  proper  earth  fixed/inertial  ref¬ 
erence  frame  transformations  and  the  correct  atmospheric  drag  acceleration  it  is 
also  necessary  to  change  the  earth  rotation  vector's  sign,  i.e.  use  -w  instead 
of  u>,  and  the  sign  of  the  drag  coefficient,  i.e.  use  -D  instead  of  D.  As  the 
time  is  stepped  backward,  the  time  for  the  nth  backward  step  given  by 

t  =  t0  -  nk  ,  n  =  1,  2,  ...,  (44) 

where  h  is  the  integration  step  size,  should  be  used  to  compute  the  associated 
luni-solar  gravitational  accelerations,  as  well  as  the  earth-fixed/inertial 
reference  frame  transformations  (see  SECTION  2).  These  modifications  apply 
to  the  variational  equations  as  well. 

As  a  result  of  the  changes  to  the  initial  conditions  for  the  equations  of 
motion,  new  variational  equation  initial  conditions  must  be  used  which  are 
consistent  with  those  for  the  equations  of  motion.  This  can  be  seen  more 
clearly  by  first  noting  the  convolutions  introduced  into  the  Keplerian  element 
set  by  making  a  velocity  sign  reversal.  Consider  Figure  9-3  which  depicts  the 
orbit  configuration  and  associated  Keplerian  elements  before  and  after  the  vel¬ 
ocity  sign  reversal.  The  following  relationships  between  the  elements  for  the 
two  cases  can  be  obtained  by  inspection: 

a' 
e' 
i* 

w' 


=  e 

-  *  -  u 

=  TT  -  to, 

=  TT  +  ft. 
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9.3.3  BACKWARD  INTEGRATION  (con't) 

To  obtain  a  new  set  of  initial  conditions  for  the  variational  equations  that 
are  consistent  with  the  equation  of  motion  initial  conditions  and  still  ref¬ 
erenced  to  the  unprimed  element  set  at  trajectory  epoch,  one  must  use  the 
following  transformation: 

!  (V  =  ♦'(to) 

to  to 


where  $  (t  )  and  $'(t  )  are  given  by  equation  (62)  of  SECTION  6,  with  the 
~.o  ~t  0 

co  o 

prime  indicating  that  the  primed  element  set  is  used  in  that  equation.  The 
K  matrix  is  given  by 


K 


I 


1  0  0  0  0 

0  10  0  0 


0  0-10 


0 


0 

0 

0 


0-10 

0  0-1 

0  0  0 


The  primed  matrix  in  equation  (48)  is  computed  using  r(tQ),  -r(t0),  and  the 
primed  element  set  in  equations  (43)  -  (54)  of  SECTION  6.  When  coordinates 


9.3.3  BACKWARD  INTEGRATION  (con't) 


9.3.4  THRJST  INTERVAL  INTEGRATION 

PULSAR  does  not  consider  thrusting  in  its  force  model  or  variational  equa¬ 
tions  and  thus  does  not  numerically  integrate  through  an  orbit  adjust  interval 
using  a  thrust  acceleration  model.  It  instead  uses  a  time  tpA  near  the  orbit 
adjust  end  time  Uq^),  where  t^  >  tQ^,  to  partition  the  trajectory  into  two 
segments.  The  first  segment  is  integrated  in  the  normal  fashion  using  the  ore- 

l 

orbit  adjust  state  vector  at  the  trajectory  start  epoch  t0  until  time  tQ^  is 
reached.  At  this  time  a  predicted  post-orbit  adjust  state  vector  is  used  as 
the  initial  state  vector  to  generate  the  second  trajectory  segment  using  the 
Cowell  integration  algorithm. 


It  should  be  mentioned  that  backward  integration  from  time  tQ/^  is  not  required, 
since  no  point  or  cross  pass  editing  will  be  done  for  some  time  period  around 
the  orbit  adjust  interval  (see  subsection  8.1.2).  Also  one  should  note  that 
cross  pass  editing  is  done  independently  on  each  trajectory  segment  and  that 
state  vector  improvements  cannot  be  propagated  across  a  segment  boundary. 


The  logical  flow  of  the  integration  path  controller  is  shown  in  Figure  9-4. 
This  controller  determines  whether  an  orbit  adjust  has  occurred  during  the 
editing  span  and  directs  the  integrator  through  the  associated  required  logic 
paths. 
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9.4  CHOLESKY’S  METHOD 

Since  6  is  real  and  symmetric,  a  Cholesky  decomposition  is  used  in  PULSAR 
to  solve  the  pass  normal  equations  (see  equations  (53)  -  (55)  of  SECTION  7). 
The  method  has  been  modified  to  identify  and  solve  for  only  those  Darameters 
which  are  well  determined. 


In  general  the  method  involves  solving  the  equation 


/»•  ' 


““t 


CA 


Aet 
~  lCA 
\  / 


~s^W 


(51) 


where  S  is  a  triangular  matrix  which  satisfies  the  following  relationships: 


S'S  =  B(tCft) 


(52) 


and 


(53) 


Note  that  in  the  equations  above  the  matrices 
S(tCA)  and  e(tCA)  are  given  by  equations  (53) 
with  the  subscripts  1 b '  and  ' e '  interchanged, 
parameter  corrections  are  given  by: 


have  been  rearranged  so  that 
-  (54)  of  SECTION  7,  respectively. 
Once  S  is  determined,  then  the 


5"  Ss<‘ca) 


(54) 
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9.4  CHOLESKY'S  METHOD  (con ‘ t) 

In  order  to  identify  those  parameters  which  are  well  determined,  the  quantity 
d  is  used  to  perform  a  relative  size  test,  where 


0  ••  for  bias  parameters, 

~  'Ot 


6.--  for  position  parameters. 


components 


Is,. 


for  velocity  parameters. 


components 


This  test  assumes  that  all  parameters  are  well  determined  unless  the  condition 


<  10‘ 


is  satisfied  for  the  YLn  parameter,  where 


-  I 


and  the  summation  vanishes  for  l  =  1.  Then  the  i  parameter  is  not  well 
determined  and  is  suppressed  as  a  solution  parameter. 

When  equation  (56)  is  satisfied,  the  parameter  is  suppressed  by  setting 
all  elements  associated  with  it  in  the  S  matrix  to  zero,  i.e. 
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CHOLESKY'S  METHOD  (con't) 


S^-  =  0  (j  =  1,  2,  •••,  Np  +  1), 


•e  Np  is  the  number  of  original  solution  parameters.  The  number  of  para- 
;rs  to  solve  for,  Ns,  is  also  reduced  by  one,  i.e. 


elements  of  the  S  and  e$  matrices  for  well  determined  parameters  are 
DUted  as  follows: 


i- 1  -!s 

I 

fe=i 


if 

s . .  ’  X  Iki  Ikj  *  ^ 


Si  -  Y. 


l  ~sk 


l- 1 

s..  -  y 

k=\ 


re,  as  before,  the  summations  vanish  for  i  =  1 .  The  inverse  of  the  S  matrix 
computed  from: 


y-i 

u/y)  -  y 


Ikl  Ikj 


\1 

--  -  .  .1 
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9.4  CHOLESKY'S  METHOD  (con't) 


where 


S’1  for  l  >_  m 
''Zm 


for  Z  <  m 


The  t^  parameter  corrections  A are  then  computed  from: 


I 


(65) 


9.5  DATA  COMPACTION 

In  order  that  any  required  post-PULSAR  processing  be  performed  more  efficiently, 
the  PULSAR  data  editor  can  compact  its  edited  and  untagged  range  difference 
tracking  data.  This  process  combines  successive  good  data  points  produced  by 
the  PULSAR  editing  processes  into  a  smaller  and  more  compact  output  data  set. 

The  data  compaction  procedure  merely  involves  adding  together  successive  un¬ 
tagged  (and  valid)  data  points  within  an  untagged  and  edited  pass  until  a  time 
span  At  from  the  time  tag  of  the  first  point  in  the  sum  is  reached,  or  until 
the  last  good  point  within  Atc  has  been  added.  This  process  is  repeated  using 
the  next  set  of  data  points  within  a  At  time  span  from  the  new  first  point 
until  the  end  of  the  pass  is  reached.  This  procedure  is  restarted  and  repeated 
on  a  pass-by-pass  basis  until  all  good  passes  have  been  processed. 
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9.5  DATA  COMPACTION  (con't) 


Also  included  with  the  compacted  data  points  are  the  associated  compacted 
observation  sigmas,  oscillator  bias  corrections,  tropospheric  corrections, 

and  antenna  offset  corrections.  To  see  how  these  compacted  corrections  are 
computed,  consider  the  compacted  data  point  P  given  by 

P  =  Z  Ct  =  ^  V[pj’  tl]  '  V{pj'  + 

■L  L 

where  V(Pj-,  t)  is  defined  by  equation  (17)  of  SECTION  7  and  the  summation 
is  carried  out  over  the  good  data  points  within  the  Atc  time  span.  The  last 
equation  can  be  rewritten  as 


The  second  and  fourth  terms  in  this  equation  are  the  compacted  oscillator 
bias  and  tropospheric  corrections  that  are  provided  with  the  compacted  data 
points.  Although  not  required  for  post-PULSAR  processing,  it  should  be  noted 
that  compacted  antenna  offset  and  general  relativity  corrections  are  obtained 
from  a  Taylor  expansion  of  the  first  term  in  the  last  equation.  Using  Equa¬ 
tion  (18)  of  Section  7,  one  may  write 
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9.5  DATA  COMPACTION  (con't) 


I  Pj  I  “  +  ArAy  *  ly  +  APQRy 


where 


X;  =  r.  -  r  . 


Thus  the  first  term  in  equation  (67)  becomes 


-*■  -*■ 
P  •  ~  P  • 


l-ll  !  =  £«,-  W  +  X 


i’ArR.  -tc-1  •  ArA .  > 


X!  (ApGR^  '  ApGRx_^  • 


The  second  and  third  sums  in  the  last  equation  are  the  compacted  -ntenna 
offset  correction  and  compacted  general  relativity  correction,  respectively. 
The  compacted  observation  sigma  is  given  by 


cc  .  <AT)-« 


where  AT  is  the  range  difference  time  interval  spanned  by  the  compaction 
summation.  The  flow  of  the  compaction  process  is  presented  in  Figure  9-5. 
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